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I. INTRODUCTION 


The original objective for this project was to demonstrate a new algorithm for synthesizing multirate sampled- 
data control laws by application to a representative aircraft control problem. That algorithm, developed in 
connection with another research effort supervised by the Principal Investigator and based on a finite-time 
quadratic cost function, eventually proved unsuitable for the aircraft control problem. To complete this project 
we therefore developed a new multirate control law synthesis algorithm, based on an infinite- time quadratic 
cost function, along with a new method for analyzing the robustness of multirate systems, and applied both to 
the aircraft control problem 

The following is a complete list of the contributions of this project: 

1 . A new generalized multirate sampled-data control law structure (GMCLS) was introduced. Features of 
this structure include an arbitrary dynamic order and structure for the processor dynamics, and sampling 
rates for all sensors, update rates for all processor states, and update rates for all actuators that can be 
selected independently, (discussed in Section II) 

2. A new infinite-time-based parameter optimization multirate sampled-data control law synthesis method 
and solution algorithm were developed, (discussed in Section HI) 

3 . A new singular-value-based method for determining gain and phase margins for multirate systems was 
developed, (discussed in Section IV) 

4. The fmite-time-based parameter optimization multirate sampled-data control law synthesis algorithm 
originally intended to be applied to the aircraft problem in this project, was instead demonstrated by 
application to a simpler problem involving the control of the tip position of a two-link robot arm. 
(discussed in Sections in and V) 

5 . The GMCLS, the new infinite-dme-based parameter optimization multirate control law synthesis method 
and solution algorithm, and the new singular-value based method for determining gain and phase 
margins were all demonstrated by application to the aircraft control problem originally proposed for this 
project (discussed in Section VI) 

These five contributions are discussed in order in the following sections of this report. The first three sections 
are in a summary form only and the reader is referred, for details, to preprints of journal papers in the 
appendixes. The next two sections present applications of the parameter optimization techniques. The final 
two sections present our conclusions and suggest topics for future research. 
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II. The Generalized Multirate Sampled-Data Control Law Structure 

A key point often ignored by the developers of multirate sampled-data control law synthesis methods is that, in 
order for any such method to be practically useful, it must provide the control law designer with the flexibility 
to independently choose the sampling rate for every sensor, the update rate for every processor state, and the 
update rate for every actuator. Such flexibility is frequently essential for efficient utilization of real-time 
control hardware, and for systems that include distributed processing and/or utilize sensors that provide only 
discrete-time signals at fixed sampling rates [1], In this section we present a general-purpose, multirate 
sampled-data control law structure (GMCLS) that provides that flexibility. 

To understand the GMCLS, it is necessary to establish a certain notation regarding the scheduling of sampling 
and update activities for a multirate system. Figure 1 shows an example of the time lines for the sampling and 
update activities of a multirate system. We define the shortest time period (STP) as the greatest common 
divisor of all of the sampling, update and delay periods; and we define the basic time period (BTP) as the least 
common multiple of all of the sampling, update and delay periods. We reserve the symbol T to represent the 
STP, and the symbol P to represent the (integer) number of STP’s per BTP. Finally, we frequendy make use 
of a doubly-indexed independent (time) variable, so that, for example, x(m,n) represents * at the start of the 
(n+l)th STP of the (m+l)th BTP, for m=0,l, ... and n=0,l, ... ,P-1. 

A block diagram of the GMCLS is shown in Figure 2. y represents the incoming, noise-free, continuous- 
time sensor signal; v is the discrete-time sensor noise signal; and u is the continuous-time control signal. The 
sampling period of the one sampler is the STP of the complete system’s sampling/update schedule. The delay 
blocks are one-STP delays; and the ZOH block represents a zero-order hold. 


Tune Lines for Sampling/Update Activities: 

1 1 

T 251 

1 Time (Seconds) 

t i i i i i i 

0 4T ST ITT 1ST 2CT 24T 

T 1 U T ? I 1 T 

0 3T ST 9T 12T 1ST 1ST 2 IT 24T 

H l-a* .1 

| BTP 1 


Figure 1 Example Multirale Sampling/Update Schedule 


4 



Figure 2 Generalized Multirate Sampled-Data Control Law Structure 


A key feature of the GMCLS is its use of the switching matrices, S^(n), and S M (n), for n= 0,1, ... ,/M, 
to represent the variations in the sensor sampling, processor state update, and control update activities, 
respectively. We define a switching matrix as a binary, diagonal matrix. S y (n) is the switching matrix that 
describes the sensor sampling activities at the start of the (n+l)/A STP (of every BTP). If the ith diagonal 
element of S y (n) is 1, then the ith sensor’s signal is sampled at the start of the (/i+l )lh STP of every BTP, and 
the sampled value, with the sensor noise v added, is immediately stored as the ith element of y. If the ith 
diagonal element of S y (n ) is 0, then the same element of y is simply held at those instants. The update activities 
for the processor state vector z and for the actuator hold state vector 5, in Figure 1, are similarly represented 
by the switching matrices S z (n\ and S u (n\ respectively, for n= 0,1, ... ,/M. 

For a detailed discussion of the GMCLS see [1]. The key points are: 

1 . The switching matrices Sz(n ), and S u (n) are completely determined by the system’s sampling and 
update activities schedule. 

2. The only unknowns are the processor matrices A z (n), B z {n\ C z (n), and D z {n) 

3 . The dynamic order of the processor dynamics (i.e., the dimension of z) is arbitrary. 

For design purposes, the implications of these points are the following: 

1 . The GMCLS provides complete flexibility with regard to the selection of sampling rates for all sensors, 
update rates for all processor states, and update rates for all actuators. The single constraint is that the 
ratio of all sampling, update and delay rates must be rational, so that the complete sampling/update 
schedule is periodic. 

2. The GMCLS provides complete flexibility with regard to the dynamic order and structure of the control 
law; i.e., the input-output dynamics of virtually any multirate sampled-data control law of practical 
interest can be realized with the GMCLS. 

3. Apart from the (significant) problem of choosing sampling and update rates, the GMCLS reduces the 
control law synthesis problem to one of determining the processor matrices A z (n), B z (n ), C z (n), and 
D z (n), for n=0,l, ... ,/M. 
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For the purpose of numerically determining A z (n), B z (n), C z (n), and D z (n) it is convenient to represent the 
GMCLS in the following state model form (see [1] for details): 


where 


c(m,n+ 1) = A c (n)c(m,n) + B c (n)y(msi) 
u(m,n) = C c (n)c(mji) + D c {n)y(mji) 

c(m,n) = [ z(m,n ) y(m,n) u(m,n)] T 

[I-SJ.n)]+S£n)A z (n) S z (n)B z (n)[I-S y (n)] 0 
0 I-S y (n) 0 

S u (n)C z (n) S u (n)D z (n)[IS y (n)] I-S u (n) J 


A c (n) = 


( 1 ) 

( 2 ) 

(3) 

(4) 


B c (n) = 


S z (n)B z (n)S y (n) 

S y (n) 

LS M (n)D z (rt)S y (n)J 


(5) 

( 6 ) 
(7) 


C c (n) = [S u (n)C z (n) S u (n)D z (n)[I-S y (n)} I-S u (n)] 

D c (n) = [5 u (n) D z {n) S y {n)]y{mji) 
with «(t) = u(mji) for all t on [(mP + n)T, ( mP + n + 1)70. 

The compensator parameters, A z (n), B z (n), C z (n), and D z (n), can be separated from the sampling schedule, 
S u (n), S y (n), S z (n), in an output-feedback representation of the GMCLS. Assuming a discretized model of 
the plant dynamics of the form 


p(mji+ 1) = A p p(m,n) + B p u(m,n) + E p v(m/i) 
y{mj i) = C p p(mji) + F p w(mji) 


( 8 ) 

(9) 


where v and w represent represent process and measurement noise, respectively, we can rewrite the closed 
loop system in the output feedback form 


( 10 ) 


- p(m,n+l) ' 


A p 0 1 

- p(mji) -j r B p 0 ir u(m,n) ' 


1 

o 

t 

- v(m,n) - 

_ c(m,n+ 1) . 


. 0 0 . 

. c(mji) J L 0 / J L c(mj t+1) - 

T 

0 0 J 

_ w(mjt) . 


y(m,n) 

L c(m,n) J 


r c p o i 

0 / 


p(mji) 

c{mji) 


i 0 f p ir 

+ 

J L o o J L 


v(mj ») 


■ u(mji) ' 


- D^n) C c (n) - 

- p(jnji) - 

. c(m,n+ 1) . 


- B c (n) A c (n ) _ 

. c(mji) . 


(ID 


( 12 ) 
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(13) 


Now the compensator 


matrices can be factored as follows. 


r D M) C c (ri) 
B c (n) A c (n ) J 


= Sl(n) 


DM) 
L BM) 


CM) - 
AM) - 


Si(n) + SM) 


where S\(n), S 2 (n), and S 3 (n) are functions of S u (n)> S y (n), and S z (n). 

Equation (1 3) is important because it allows us to separate the unknown compensator parameters A z (n), B z (n ), 
C z (/i), and D z (n) from the known sampling schedule. 

In the following section we will introduce two synthesis algorithms that can be used to determine the optimum 
compensator parameters A z (n), B z (n), C z (n), and D z (n). 


III. Parameter Optimization Control Law Synthesis Methods 

There are five well-recognized techniques for synthesizing multirate control laws: successive loop closures, 
pole placement, singular-perturbation-based methods, LQG Optimal methods, and parameter optimization 
methods. 

The advantages of successive loop closures are that it is easy to use, that it can (conceivably) be used to 
synthesize control laws of arbitrary dynamic order and structure, and that it is particularly effective in 
applications where the control loops are not strongly dynamically coupled. Its disadvantage is that its one- 
loop-at-a-time approach cannot fully account for all dynamic coupling between control loops. 

The problem with pole placement is determining where the closed-loop poles should be placed. It is a 
particularly difficult problem in the multirate (as compared to the single-rate) case because the STP-to-STP 
dynamics of multirate systems are periodically time-varying [2]. Only the BTP-to-BTP dynamics of multirate 
systems are time-invariant, and it is the poles of those dynamics that are assigned by pole-placement. In 
applications, determining desirable BTP-to-BTP closed-loop poles for a typical multirate system is difficult 
because the BTP of its sampling/update schedule will typically be longer than many of its desired closed-loop 
characteristic times. 

Singular-perturbation-based control law synthesis methods amount to successive loop closures prefaced with a 
coordinate transformation to separate the full control law synthesis problem into two or more dynamically 
decoupled control law synthesis problems of different time scales. A complete decoupling requires changes in 
not just the state coordinates, but in the input and output coordinates as well. Such a decoupling is not 
possible in the multirate case because the input and output coordinates represent physical sensor and actuator 
signals destined to be sampled/updated at different rates. 

The advantages of the LQG optimal control law synthesis methods are that stabilizing control laws are 
relatively easy to obtain and that the control laws for all control loops are synthesized simultaneously, taking 
full advantage of all dynamic coupling between the control loops. The disadvantages are that the dynamic 
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order and structure of the control law is fixed, that stability robustness objectives are difficult to achieve, and 
that the resulting control laws are periodically time-varying [2]-[3]. 

We favor parameter optimization methods for control law synthesis for multirate systems because they offer 
the principal advantages of the successive loop closures and LQG optimal synthesis methods. These 
advantages are that control laws of arbitrary dynamic order and structure can be synthesized, and that control 
laws for all control loops can be synthesized simultaneously, taking full advantage of all dynamic coupling 
between control loops. The disadvantage of parameter optimization methods is that a numerical search is 
required to determine the control law parameters. 

In this section we present two parameter optimization methods for synthesizing multirate control laws. Both 
utilize the GMCLS discussed in Section II. The first is based on a Finite-time quadratic cost function while the 
second is based on an infinite-time quadratic cost function. Both methods solve the multirate compensator 
synthesis problem by using a gradient-type numerical search to find a set of compensator parameters that 
minimize a quadratic cost function. 

The multirate optimization problem is as follows. 

Given: 

1. The plant dynamics represented by 

p (0 = Ap pit) + B pu m +Bpv v (0 ( 14 ) 

y(t)-C p ~p(t) 05) 

Here p is the plant state vector, u is the control input vector, y is the noise-free measurement output 
vector, and v is the noise input vector. 

2. The complete sampling and update schedule for the compensator. This amounts to specifying S u (n), 
Sy(n), and S z (n), for n=0,l, ... , P-1. 

3 . The order for the processor dynamics (the number of elements in z in (3)). 

4. The desired structure (e.g., a diagonal structure) for the processor matrices, A z (n), B z (n), C z (n), and 
D z (n), for n=0,l, ... , P-1. 

5 . The number of distinct sets of processor matrices and when they are active. The optimization algorithms 
allow A z (n), B z (n), C z (n), and D z (n ) to be periodically time varying. The designer can specify equality 
relations among the compensator matrices. For example, if a time invariant compensator is desired then 
the designer can specify that A z (0) = A z (l) = ••• = A Z (P-1), and similarly for B z ,C Z and D z . 

6. The power spectral density V of the process noise v (in (8)). 
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7. The covariance Win), for n=0,l, ... , P-1, of the sensor noise w (in Fig. 1). w is assumed to be a 
periodically stationary, gaussian, purely random sequence, with period equal to the BTP of the 
sampling/update schedule. 


8 . The time tf and non-negative definite weighting matrices Q and R for the performance index 


%) = 


f 

_L f 

'pit) 

T 

Q o' 

2 t f J 
l f o 

_«(0. 


.OR. 


'pit) 


dt > 

J uit )_ 

> 


(16) 


where E is the expected value operator. 

In the finite time optimization problem /ymust be a multiple of the BTP of the sampling/update schedule. 
In the infinite time optimization problem fy-— »°° and J infinite-time ~ t lim J(t f ) 


Find: 

A set of processor matrices, A z {n), B z (n), C z (n), and D z (n), for n=0,l, ... ,P- 1, such that the performance 
index 

J ss = 

tf—t <» 


is minimized. 

This optimization problem can be solved using either the finite-time cost function or the infinite-time cost 
function. 

Solution Method Using the Finite-Time Cost Function, 

The finite-time optimization algorithm was developed in connection with another research effort supervised by 
the Principal Investigator. This method synthesizes the multirate compensator that minimizes /(*/) for a finite 
A detailed discussion of this method can be found in [1]. A summary of the solution procedure follows. 

1 . Determine closed-form expressions for the performance index J{tfi> and for its gradients with respect to 
the elements of the processor matrices A z (n\ B z (n ), C z {n\ and D z (ri) , for n- 0,1, ... */M. 

2. Use a gradient-type numerical optimization algorithm to determine a set of processor matrices, A z {n ), 
5 z (n), C z (n), and D z {n ) , for n=0,l, ... ^-1, that minimizes 

3 . Obtain a steady-state solution by re-optimizing for larger and larger ryuntil tf gets to be large compared to 
all of the closed-loop system’s characteristic times. 

The advantage of this method is that with tf finite, the cost function J(tf) remains finite even if the compensator 
is destabilizing. The designer does not need to find a stabilizing compensator to start the optimization process 
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as long as If is small enough that Jiff) does not exceed the numerical limits of the computer performing the 
optimization. 

The disadvantage of this method is that the closed-form expressions that have been developed thus far for the 
performance index J(tj) and for its gradients with respect to the elements of the processor matrices are very 
complex and computationally intensive. In addition, we encountered difficulties when applying this method to 
the aircraft control problem because our solution algorithm lacked provisions for automatic scaling of the 
control law parameters (i.e., the independent variables) during the numerical search. The sheer complexity of 
the finite-time performance index and gradient expressions prevented us from adding the automatic scaling 
provisions that would have allowed us to apply this method to the aircraft control problem. 

Solution Method Using In finite-Time Cost Function 

Instead of modifying our existing finite-time-based algorithm to alleviate the scaling problem discussed in the 
previous paragraph, we chose to develop a new infinite-time-based multirate sampled-data control law 
synthesis method, based on corresponding developments for single-rate systems by Mukhopadhyay [4], for 
which much simpler performance index and gradient expressions are easy to derive. For a complete 
description of that method, and the solution algorithm we developed to implement it see [5]. A summary of 
the solution procedure follows. 

1. Find an initial stabilizing guess for the processor matrices A z (n) f B z (n), C z (n), and D z in) f for 
n= 0,1, ... ,/M. The finite-time solution algorithm requires an initial stabilizing compensator because 
J ss is infinite when the closed loop system is unstable. From our experience, many multirate problems 
can be stabilized using successive loop closures. The aircraft problem was open loop stable, and so 
determining a stabilizing compensator was trivial. 

2. Determine the necessary conditions (given in [5]) for the processor matrices, A z (n), B z {n) y C z (n), and 
D z (n), for *=0,1, ... ,/M to minimize J sr These are represented by three sets of coupled matrix 
equations. Two sets are Lyapunov equations, one governs the steady state covariance of the plant and 
control states, and the other governs a Lagrange multiplier. The third represents the gradient of J ss with 
respect to the compensator parameters. 

3 . Use a gradient- type numerical search to solve the necessary conditions and determine a set of processor 
matrices, A z (n), B z (n\ C z (n), and D z (n), for *=0,1, ... ,/M, that minimizes 

The advantage of this method is that the gradient of J ss with respect to the compensator parameters is easy to 
evaluate via the necessary conditions. For a given problem, the infinite-time optimization algorithm typically 
requires fewer computations to find the optimum compensator parameters than does the finite-time 
optimization algorithm even when both algorithms are initialized with the same stabilizing compensator. 

Even though the finite-time and infinite-time based solution algorithms can determine optimum compensator 
parameters, there is no guarantee that the design will be robust. In the following section we present a method 
for analyzing the robustness of a multirate control system. 
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IV. Gain and Phase Margins for multirate Systems Using Singular-Values 

There are many established methods for synthesizing multirate compensators, see Section IE, but surprisingly 
few methods for analyzing the robustness of these systems. Current robustness analysis methods rely 
principally on the transfer function of the system. A multirate transfer function, in the traditional sense, does 
not exist, because multirate systems are periodically time varying. Without modification, established single- 
rate analysis methods cannot be applied directly to multirate systems. 

As part of this project, we developed an approach for extending the nyquist criterion and singular value 
analysis to multirate and periodically time varying systems. For a detailed discussion of this approach, 
including application of structured singular value robustness analysis to multirate systems, see [6]. In this 
section we present a summary of the important ideas from that paper used to calculate gain and phase margins 
of multirate systems using singular values. 

As we saw in Section II, a multirate compensator can be modeled as a linear periodically time varying system 
(l)-(2). Equations (l)-(2) from Section II can be written as 

c(m,n+l ) = A c (n)c(mjt) + B c (n)y{mji) (17) 

u{m,n) = C c (n)c(mji) + D c (n)y(mjn) ( 1 8) 


This system (17)-(18) can then be transformed to an equivalent single-rate system (ESRS) by repeated 
application of (17)-(18) over the BTP [7]. The ESRS has the form: 

c(m+l ,0) = A e c(m, 0) + B e y{m, 0) ( 19) 

H(m,0) = C e c(m,0) + Djy(m, 0) (20) 


where y(m, 0) = 


?(/n,0) 

y(«.l) 

y(mf-n 


u(m, 0) 
u(m, 1) 

u(mS- 1) 


The transfer function for the ESRS is 

u(z p ) = G P (z p )y{z p ) (22) 

where Gp(z p ) = C e (/z^ - A e )^B e + D e (23) 

For a detailed discussion of the ESRS, see [6]. The key points are: 

1 . The ESRS is a time invariant single-rate system with a sampling period of one BTP and the unique 
property that the inputs are time correlated and the outputs are time correlated. 
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2 . In general G P (z p ) has a very complicated form, but it can be shown that if the system is time invariant 
with G(z) equal to a constant, then G P (z p ) will also be constant and block diagonal with G(z) on the 
diagonal. 

3. The ESRS allows us to manipulate time invariant and periodically time varying systems (e.g. multirate) 
as if they were both time invariant. The state space or transfer functions descriptions can be used to 
calculate input-output relations for systems in series or in a feedback loop just as in classical control [8]. 
For example, to calculate the ESRS of a multirate compensator in series with a time invariant plant, we 
would calculate the ESRS of the plant and compensator individually and then combine them using block 
diagram arithmetic. 

4. Kono [9] has shown that if the ESRS is stable then the muitirate system from which it was derived will 
be stable. 

5. Single-rate robustness analysis techniques can be applied to the ESRS as long as the results are 
interpreted in light of the fact that some of its inputs and outputs are time correlated. 

Generalized gain and phase margins for the ESRS (and equivalendy the multirate system) can be calculated 
using singular value analysis. If we assume a plant uncertainty of the form 

G<, 2 ) Actual = G(z) Nomimd kei Q (24) 

then the ESRS plant uncertainty has the form 

G p (z p ) Actual - G p (z p )^ om i na i(.lcei & ) P (25) 

(ke/ Q ) P = diag [kei Q , ke ^, . . . , ke&] with P blocks 
[Recall that if H(z) is constant, H P (z p ) is block diagonal with H(z) on the diagonal.] 

The multirate system is guaranteed to remain stable whenever 

^((fce/ 6 )* 1 -1) < <?(/ + G P (z p )) on the nyquist contour (26) 

Traditional gain margins can be obtained by setting 0 = 0 and solving (26) for k. Phase margins can be found 
by setting k = 0 and solving (26) for 0. 

As with most singular value robustness analysis methods, the k and 0 found using (26) are conservative. If, 
however, kefi is diagonal, the conservativeness associated with (26) can be reduced by diagonally scaling 
G P (z p ). We used Osborne’s method of preconditioning matrices to increase the lower bound for G P {z p ) and 
thus to improve our estimate of the gain and phase margins. 
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V. Application of the Finite-Time-Based Parameter Optimization algorithm to 

a Two Link Robot arm Control problem 

The original proposal for this project called for the finite-time-based parameter optimization multirate sampled- 
data control law synthesis method of Section IE to be applied to an aircraft control system design problem. 
That method and a solution algorithm to implement it had been previously developed as part of another 
research effort supervised by the Principal Investigator. Due to the solution algorithm’s lack of adequate 
provisions for automatic scaling of the control law parameters (i.e., the independent variables) during the 
numerical search, we were not able to apply it successfully to the aircraft control problem. We maintain, 
however, that the problems we encountered with it were a consequence of problems with the solution 
algorithm and are not necessarily indicative of problems with the synthesis method. 

In this section we therefore present an application of the finite-time-based multirate sampled-data control law 
synthesis method to a two-link robot arm control problem. The robot arm application demonstrates the utility 
of the method without being so poorly conditioned that automatic scaling of the control law parameters was 
required during the numerical search. 

The two-link robot arm system we dealt with is shown in Figure 3. The first link is long and massive, for 
large-scale slewing motions. The second is short and lightweight so high-bandwidth control of the tip position 
can be achieved with a relatively small motor at the second joint. The pin joint, rotational spring, and 
rotational damper at the midpoint of the first link model flexibility in that link. The control inputs are the motor 
torques, T\ and T 2 . The measured outputs are the joint angle 0 and the tip position 8. The spring constant (k) 
and damping coefficient ( b ) values (in Fig. 3) yield an open-loop vibration mode with a 10 Hz natural 
frequency and 1% damping. 



L\ 0.5 kg 0.5 m 

L 2 0.5 kg 0.5 m k = 37.33 N/rad 

L 3 0.04 kg 0.2 m b = 0.01 2 N s/m 

The natural frequency of the vibration mode is 10 hz. 
Inputs: Torques T i and T 2 
Outputs: 9 and 5 

Figure 3 Two-Link Robot Arm System 
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We used the finite-time-based multirate control law synthesis method of Section HI to synthesize multirate 
sampled-data control laws for this system. Our performance objective was high-bandwidth control of the tip 
position 8, and it is intuitively clear that this can best be accomplished, given a fixed real-time computation 
capability, by trading low-bandwidth control at Xj for high-bandwidth control at T^. Thus, for an 8-to-l 
control bandwidth ratio, we chose the sampling/update rate for 8 and to be 8-times faster than that for 8 and 
7\. For comparison purposes, we designed also corresponding analog and single-rate sample-data control 
laws. 

TLA Control Laws 

For the TLA system, the tip position (8) responses to a commanded step change in the tip position obtained 
with the analog, single-rate and multirate control laws we synthesized are shown in Figure 4. See [1] for 
additional results and details. A summary description of those designs follows. 

LQR Analog Design The LQR Analog response was obtained with an analog LQR (full state feedback) 
control law that is optimal with respect to a quadratic performance index that yields 0.7071 damping 
= 0.7071) and an 8-to-l ratio of characteristic frequencies (a)„ 2 /o ) nl = 8) for the two closed-loop 

modes. 

Third-Order Analog Successive Loop Closures Design The Third-Order Analog Successive Loop Closures 
response was obtained with a successive loop closures control law that consisted of a single lead network from 
0 to Tj, and two identical cascaded lead networks from 8 to T 2 . The gains, and zero and pole locations were 
chosen to yield dominant closed-loop poles coincident with those obtained with the LQR Analog control law. 

Third-Order Multirate Tustin Design The Third-Order Multirate Tustin response was obtained with a 
multirate sampled-data control law obtained by discretizing the lead compensators of the Third-Order Analog 
Successive Loop Closures design using Tustin’s method [10]. The 9-10-7] control-loop sampling/update rate 
is a factor of 8 times the characteristic frequency of the lower-frequency closed-loop mode from the Third- 
Order Analog Successive Loop Closures design; and the 8-to-72 sampling/update rate is the same multiple of 
the characteristic frequency of the higher-frequency closed-loop mode from the Third-Order Analog 
Successive Loop Closures design. 

Optimized Third-Order Multirate Tustin Design The Optimized Third-Order Multirate Tustin response was 
obtained with a control law synthesized by the finite-time-based multirate sampled-data control law synthesis 
method of Section HL This control law is the Third-Order Multirate Tustin control law, but with its gains and 
its pole and zero locations optimized to minimize the same performance index as is minimized by the LQR 
Analog control law. 

Analog Third-Order Design The Analog Third-Order response was obtained with a third-order, generally- 
structured, analog control law synthesized using Ly’s Sandy algorithm [11]-[12] to minimize the same 
performance index as is minimized by the LQR Analog control law. 
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Multirate First-Order, Second-Order & Third-Order Designs The Mukirate First-Order, Second-Order, and 
Third-Order responses were obtained with muitirate, generally-structured, sampled-data control laws 
synthesized by the finite-time-based multirate sampled-data control law synthesis method of Section III to 
minimize the same performance index as is minimized by the LQR Analog control law. The sensor sampling 
and actuator update rates are the same as in the Third-Order Multirate Tustin control law. In the First-Order 
case, the update rate for the one processor state is the same as the faster sensor-sampling/actuator-update rate. 
In the Second-Order case, one processor state is updated at the faster rate and the other at the slower rate. In 
the Third-Order case, two processor states are updated at the faster rate and one is updated at the slower rate. 

Single-Rate Third-Order Design Finally, the Single-Rate Third-Order control law response was obtained 
with a single-rate, generally-structured, sampled-data control law synthesized by the finite-time-based multirate 
sampled-data control law synthesis method to minimize the same performance index as is minimized by the 
LQR Analog control law. Its single sampling/update rate was chosen to require the same average number of 
computations per unit time for real-time operation as is required for real-time operation of the Multirate Third- 
Order control law. 

Conclusions 

The TLA results in Figure 4 demonstrate some of the benefits of multirate control. For example, the tip 
position overshoot (5) with the multirate compensator is much less than with its equivalent single-rate 
counterpart But more importantly, the results demonstrate that the finite-time-based multirate sampled-data 
control law synthesis method can be used to synthesize multirate control laws of arbitrary structure and 
dynamic order, with arbitrarily selected sampling rates for all sensors, and update rates for all processors states 
and actuators. The third-order multirate compensator, for example, uses two different update rates for the 
processor states, inputs and outputs, and a general compensator structure with full coupling between inputs, 
outputs and processor states of different rates. 
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VI. Application of the Infinite-Time-Based Parameter Optimization algorithm 
to a Yaw Damper and Modal Suppression System for a Commercial aircraft 


A practical application of multirate control can be found in aircraft. The limited computational resources of 
aircraft dictate that their control systems must function efficiently. Multirate control allows the designer to 
efficiently allocate these resources by trading slow sampling and update rates in control loops associated with 
low-bandwidth control functions for fast sampling and update rates in control loops associated with high- 
bandwidth control functions. In this section we consider a particular application of multirate control, a 
combination yaw-damper and modal suppression system for a commercial aircraft. 


In the interest of weight reduction for fuel efficiently, aircraft are being constructed with less structural rigidity. 
Structural vibration modes can be excited in such aircraft by wind gusts or by movements of control surfaces. 
These vibrations affect not only the structural integrity of the fuselage but also passenger ride quality. In the 
lateral direction, such vibrations are often induced by rudder activity associated with the yaw-damper. A 
“modal suppression system” can be added to the yaw-damper loop to suppress these vibrations. The modal 
suppression system would traditionally be designed by successive loop closures. 


In this section we describe the design of a multirate combination yaw -damper and modal suppression system 
for a commercial aircraft using the infinite-time-based multirate compensator synthesis algorithm and 
robustness analysis technique discussed in Sections III and IV. For comparison purposes we also designed 
corresponding analog and single-rate sample-data systems. 

The goal for each compensator design was to increase the damping of the dutch-roll mode to 0.6, and to 
decrease the covariance of lateral accelerations at the nose and aft of the airplane, particularly those components 
associated with low frequency flexible modes. The performances of the compensators were compared by 
comparing the closed loop dutch-roll damping, the covariances of lateral accelerations at the nose and aft of the 
aircraft due to a unit covariance gaussian white noise disturbance, and the PSD plots of the lateral accelerations 
at the nose and aft of the aircraft for either a white noise disturbance (analog designs) or a gust pulse 
disturbance (sampled-data designs). 

Open Loop Aircraft 


A block diagram of the airplane model is shown in Figure 5. The lateral dynamics model consists of 4 rigid 
body modes (heading, spiral, dutch roll and roll) and 1 1 flexible modes. Actuator/power control units for the 
aileron and rudder are modeled as second-order lags. 


30(35) ___ 

^S)Rudder "(j + 30)(j + 35) 


20(25) 

G ( S ) Aileron = (j + 20)(j + 25) 


(27) 


The lateral gust disturbances are filtered by a second-order Dryden gust model 


G(s) = 


17.4965 + 2.1617 
21.836s 2 + 9.3458$ + 1 


(28) 
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Figure 5 Block Diagram of the Open Loop Airplane Model 


The poles associated with the spiral and heading modes were compensated with static gain feedback before the 
yaw-damper/modal suppression systems were designed, because these modes, which laid close to the origin 
and were controllable with the rudder, created numerical difficulties for Sandy (the optimization program used 
to design the Fourth-Order Analog compensator discussed in later in this section). The spiral mode was 
compensated by feeding back roll and roll-rate to the aileron. Heading was compensated with heading to 
rudder feedback. In what follows we refer to the airplane model with spiral and heading modes compensated 
as the uncompensated airplane (no dutch-roll compensation). 

The lateral accelerations of the uncompensated airplane are measured by Nynose and Nyaft. The PSD plots of 
lateral accelerations for the uncompensated airplane are shown in Figure 6. A yaw-damper/modal suppression 
system should reduce the total area under this curve (covariance of lateral acceleration). In particular, it should 
reduce the peak at * 0.5 Hz (near the dutch-roll mode) and the peaks between 3 Hz and 6 Hz (low frequency 
flexible modes). Values of the dutch-roll damping, and the Nynose and Nyaft covariances for the 
uncompensated airplane are given in Table 1. 
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Figure 6 PSD of Nynose and Nyaft for Uncompensated Airplane with Unit 
Covariance Gaussian White Noise Lateral Disturbance 


Table 1 Results for Analog Designs with a Unit Covariance Gaussian White Noise 
Lateral Disturbance 


Design 

Dutch-Roll Damping 

1 

Nynose Cow. 

(ft 2 /sec 3 ) 

Nyaft Cov. 
(ft 2 /sec 3 ) 

Uncompensated 

0.08 

5.1 

21.8 

j Analog Yaw-Damper Only 

0.6 

5.0 

6.1 

i LQR Analog 

0.6 

2.4 

3.1 

! Fourth-Order Analog 

0.55 

2.5 

2.4 


Analog Yaw-Danroer/Moda l Suppression System Designs 

Three analog compensators were designed: a yaw-damper only system, a full state feedback yaw- 
damper/modal suppression system, and a fourth-order yaw-damper/modal suppression system. PSD plots of 
Nynose and Nyaft for the analog designs are shown in Figure 8; Nynose and Nyaft covariances are 
summarized in Table 1 . These analog designs provide a base line for comparison with the sampled-data 
designs and were used to determine appropriate values for cost weighting matrices. Following is a summary 
of these designs. 
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Analog Yaw-Damper Only Design The yaw-damper only design uses static feedback from y to 8,, using a 
gain k damper We chose k damptr such that the dutch-roll damping was 0.6 using classical root locus. While the 
peak on the PSD plot associated with the dutch-roll mode (=0.5 Hz) has been reduced significantly from the 
uncompensated case, the peak near 3 Hz has increased (see Fig. 8). This is the problem with using static gain 
feedback. As you “press” on one peak of the PSD another “pops” up due to the input coupling between the 
dutch-roll and low frequency flexible modes. 


LQR Analog Design The LQR design uses full state feedback to improve the dutch-roll damping and reduce 
the covariance oi Nynose and Nyaft. The compensator was designed to minimize the following cost function. 


J ss = lim E l 

tj — » «> 


Nynose 
L Nyaft J 


0.001 0 
. 0 0.004J 


Nynose 
L Nyafl J 


+ 1 . 68 . 


(29) 


Weighting matrices for (29) were chosen such that the covariances of Nynose and Nyaft were reduced from 
the yaw-damper only case by the same percent, and the dutch-roll mode had a damping of 0.6. Figure 8 
shows that the LQR design significandy reduces the dutch-roll peak as well as the peaks associated with the 
flexible modes. 

Fourth-Order Analog Design The Fourth-Order Analog compensator is a yaw damper/modal suppression 
system designed using Sandy [11]. A block diagram of this compensator is shown in Figure 7. This design 
mi nimise the same cost function as the LQR design (29) with the weighting on 8 r adjusted to achieve close to 

the desired 0.6 dutch-roll damping. An unexpected result is that the covariance at Nyaft for the Fourth-Order 
Analog design is actually better than for the LQR Analog design. This is a consequence of adjusting the cost 
function weighting matrices to achieve the desired the dutch-roll damping. 
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Figure 7 Block Diagram of Airplane with Fourth-Order Analog Compensator 
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Figure 8a PSD of Nynose for Aiiplane with Analog Compensators for Unit 
Covan ance Gaussian White Noise Lateral Disturbance 
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Sampled-Data Yaw-Damper/Modal Sup pression System Designs 

Three sampled-data compensators were designed: a single-rate yaw-damper only system, a fourth-order 
multirate compensator and a fourth-order single-rate compensator. Both fourth-order compensators were 
synthesized using our infinite-time-based multirate control law synthesis algorithm to minimize the same cost 
function as the LQR Analog design. 

The sampled-data compensator designs were based on a maximum sample/update rate of 50 Hz. This is 10 
times the rudder actuator roll off frequency and 8 times faster than the fastest flexible mode which contributes 
signiftcandy to the PSD of the lateral acceleration. This sample rate is close to the slowest practical sample rate 
which could be used. 

PSD plots for the sampled-data designs were generated using a gust pulse (a rectangular pulse) at the 
disturbance input, as opposed to the gaussian white noise used for the analog designs. For the analog 
designs, the PSD plots were based on transfer functions from the disturbance input to Nynose and Nyaft. 
Multirate compensators are periodically time varying so that transfer functions for them, in the traditional 
sense, do not exist. For this reason, we used the gust pulse disturbance input to generate the PSD plots for the 
sampled-data designs. 

The gust pulse input PSD has a connection to the white noise input PSD. For a time invariant continuous 
system, the PSD plots generated using either gaussian white noise or a continuous impulse input are exacdy 
the same. This is because the Fourier transform of the impulse response is the same as the bode plot, and the 
PSD of gaussian white noise is a constant. If a continuous system, given by x(l) = Ax(t) + Bu(t), is such that 

t p 

B = u Je AT Bdr ( 30 ) 

where u can be selected arbitrarily and T p is much shorter than the observation time, then a continuous impulse 
can be approximated by a pulse of duration T p and magnitude u. For the airplane problem addressed in this 
project, (30) is satisfied for 

T p = 0.02 seconds and u- 50 ft/sec. (31) 

PSD plots of Nynose and Nyaft for the sampled-data designs are shown in Figure 11. Nynose and Nyaft 
covariances for these designs are summarized in Table 2. Following is a summary of the sampled-data 
designs. 

Single-Rate Yaw-Damper Only Design The Single-Rate Yaw-Damper Only design is similar to the Analog 
Yaw-Damper design except that a sampler is used at the output y and a zero order hold is used at the input 5 r . 
Both the sampler and zero order hold operate at 50 Hz. The performance of the sampled-data yaw damper is 
very close to that of the analog damper (Figs. 8 and 1 1). 
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Multirate Fourth-Order Design The multirate compensator is shown in Figure 9. It was designed to 
minimize the same cost function as the LQR Analog design with the weighting on 8 r adjusted to achieve the 

desired dutch-roll damping. The compensator uses two sampling/update rates. The rudder is updated and the 
lateral accelerations are sampled at 50 Hz; \\f is sampled at a slower rate, 12.5 Hz, because it is composed 
primarily of the slow dutch-roll mode. 

Two of the processor states for this multirate compensator are updated at the fast rate, 50 Hz, and two are 
updated at the slow rate, 12.5 Hz. Initially a compensator was designed in which all of the processor states 
were updated at 50 Hz, but we found that there was no noticeable performance degradation if two of the 
processor states were updated at the slower rate. Slowing the update rate of these states reduces the number of 
computations required per unit time for real time implementation of the multirate compensator. 


Table 2 Results for Sampled-Data Designs with a Unit Covariance Gaussian White 
Noise Lateral Disturbance 



Figure 9 Block Diagram of Airplane with Multirate Fourth-Order Compensator 
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Single-Rate Fourth-Order Design The Single-Rate Fourth-Order compensator is shown in Figure 10. The 
sampling rate for this compensator is 28.6 Hz. That rate was chosen such that the number of multiplications 
required per unit time for its real time operation is the same for the multirate compensators. The cost function 
used to design the single-rate compensator was the same as was used to design the Multirate Fourth-Order 
compensator. 



Figure 10 Block Diagram of Airplane with Single-Rate Fourth-Order Compensator 
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Figure 11a PSD of Nynose for Airplane with Sampled-Data Compensators for 
Gust Pulse Lateral Disturbance 






Gain and Phase Margins for Sampled-Data Designs 


Gain and phase margins at the control input (8 r ) were evaluated for the Multirate Fourth-Order and Single-Rate 
Fourth-Order compensators using the robustness analysis methocs of Section IV. Table 3 summarizes the 
traditional gain and phase margins for the these compensators. Figure 12 shows the region of guaranteed 
stability for simultaneous changes in k and 8 for both compensators. 


Table 3 Traditional Gain and Phase Margins 


Design 

Gain Margin (db) [0 = 0} 

Phase Margin (Deg) [k = Odb] 

Multirate Fourth-Order 
Single-Rate Fourth-Order 

[-3.8, 7.1] 
f-3.3, 5.5] 

±32* 

±27* 



Figure 12 Stability Region for Simultaneous Gain and Phase Uncertainty for Fourth- 
Order Sampled-Data Compensators 
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Conclusions 


Figures 8 and 11 show that the yaw -damper/modal suppression systems significantly decrease the covariance 
of the lateral acceleration at the nose and aft of the airplane while attaining the desired 0.6 dutch-roll aamping. 
It should be no surprise that the analog compensators out performed the sampled-data compensators because 
the sampled-data compensators were designed using a slow sampling rate. Still, both fourth-order sampled- 
data compensators reduce the peak accelerations of Nynose and Nyaft by 175% and 50% respectively over the 
yaw-damper only systems. The performance of the Single-Rate Fourth-Order compensator is nearly as good 
as that of the multirate compensator, but, for input gain and phase uncertainty, the multirate compensator is 
more robust than the single -rate compensator. 


VII. SUMMARY AND CONCLUSIONS 

In this report we have presented a methodology for designing multirate control systems. We have introduced 
the Generalized Multirate Control Law Structure (GMCLS) which allows complete flexibility with regard to 
the dynamic order and structure of the control law, and with regard to the sampling rates for all sensors and the 
update rates for all processor states and actuators. We have presented two parameter optimization multirate 
control law synthesis algorithms, one based on an infinite-time cost function and the other based on a finite- 
time cost function, which can be used to find optimum values for the GMCLS parameters. We have presented 
a technique for determining gain and phase margins for multirate systems. Finally, we have demonstrated our 
methodology by applying it to the design of a two link robot arm control system and to the design of a 
combination yaw-damper and modal suppression system for a commercial aircraft. The application to the 
aircraft control problem, in particular, demonstrates that the methodology can be applied to design problems of 
a scale that one might expect to encounter in practice. 


VIII. Suggestions for Future Research 

The results presented here demonstrate a methodology for multirate digital control system design that is 
applicable to practical problems. Before this methodology can be routinely applied in practice, however, the 
following need to be developed: 

1 . A means for directly synthesizing robust multirate control laws. 

2. Numerical optimization algorithms incorporating auto-scaling of the independent variables and other 
features that more effectively deal with the practical difficulties of parameter optimization applied to 
multirate control law synthesis. 

With regard to direct synthesis for robustness, there are several possibilities. One would add the multiple- 
plant-condition design for robustness ideas of Ly [ 10]-[1 1]. A second would add direct nonlinear robustness 
constraints on the control law parameters during the numerical optimization. The lauer approach has been 


27 


successfully applied by Mukhopadhyay [4] to synthesize robust single-rate control laws by parameter 
optimization. 

In addition to theoretical work, a second major research effort in multirate control needs to be directed toward 
experimental research. Now that a bonafide multirate control system design methodology has been developed, 
we strongly believe that further substantive progress in the field can best be made in conjunction with bonafide 
hardware applications of that methodology in the laboratory. 
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I. INTRODUCTION 


Even in this age of fast, low-cost microprocessors there remain several important motivations 
for muitirate sampling in sampled-data control systems. One need only consider large space 
structure control problems to realize that the cost, bulk, and weight of real-time computing 
hardware continues to be an important control system design issue. Multirate sampling provides 
the opportunity to allocate sampling rates, and thus real-time computing power, more efficiently. 

In two-time-scale control problems, for example, muitirate sampling allows slow sampling in 
control loops associated with low-bandwidth control functions to be traded for fast sampling in 
those associated with high-bandwidth control functions. 

As with microprocessors, the costs of analog-to-digital and digital-to-anaiog converters are aiso 
computation-rate dependent. Multirate sampling thus provides another opportunity to reduce 
hardware costs because the computation rates required of analog-to-digital and digital-to-analog 
converters frequently depend upon their sampling rates. Multirate sampling can even be used to 
reduce the total number analog-to-digital and/or digital-to-analog converters required by a system, 
by sample-dependent scheduling of multiple conversion tasks to a lesser number of conversion 
devices. 

A third "motivation" for muitirate sampling is becoming increasingly important: sometimes 
multirate sampling is the only choice. This situation can arise when an apriori decision has been 
made to include in a system a sensor that provides a discrete-time signal at a fixed sampling rate. A 
head position control system for a computer disk drive is a good example of such a system. The 
disk head, which is suspended atop the rotating disk, includes a sensor that reads the head position 
directly from certain diametrically-spaced segments on the disk. The sensor's sampling rate is thus 
fixed by the disk's rotation speed. To increase the control bandwidth beyond that dictated by that 
sampling rate, a second, faster-rate sensor must be added. 

A key point often ignored by developers of multirate control law synthesis methods is that 
these motivations for multirate sampling dictate also certain flexibilities required to meet the needs 
of engineering practice. Specifically, multirate control law synthesis methods, to meet the needs of 
engineering practice, must allow the sampling rates for all sensors, the update rates for all processor 
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states, and the update rates for all actuators to be specified independently. The one generally 
accepted restriction, with regard to these rates, is that the ratio of all combinations of sampling and 
update rates must be rational, so that the complete sampling/update schedule will always be 
periodic. (We assume that all sampling and update events are synchronized to the same clock. 

The asynchronous case is treated elsewhere (1).) 

Time lines representing such a periodic sampling schedule are shown in Fig. 1. We define the 
Basic Time Period (BTP) of such a schedule as the least common multiple of all of its sampling and 
update periods. The BTP is the period of repetition of the sampiing/update schedule. We define 
the Shortest Time Period (STP) as the greatest common divisor of all of its sampling and update 
periods. We reserve the symbol P to represent the (integer) number of STP's per BTP, and we shall 
frequently use a double-indexing scheme for the independent variable so that, for example, x(m.n) 
represents z at start of the (n+l)th STP of the (m+l)th BTP, for m - 0,1 and n = 0, P - 1. 

There are five well-recognized methods for synthesizing multirate sampled-data control laws: 
successive loop closures, pole placement, the singular perturbation method, the LQG (linear 
quadratic Gaussian) method, and parameter optimization methods. Successive loop closures (2] is 
arguably the most important because it is the single one of the five that is widely used in industry. 
The advantages of successive loop closures are that its one-loop-at-a-time approach requires no 
new multirate synthesis techniques, and that the sampling/update rate for each control loop can be 
specified independently. The problem with successive loop closures is that its one-loop-at-a-time 
approach cannot fully account for all dynamic coupling between control loops. 

Pole placement [3,4,5,6,71 for multirate systems has received considerable recent attention in the 
wake of reports on the capacity for periodically time-varying output feedback controllers to place 
closed-loop poles. In Ref. 3, for example, it is shown that given any controllable and observable 
continuous-time plant with m inputs, it is always possible to construct a periodically time-varying, 
pure-gain, output feedback control law that places the closed-loop poles arbitrarily, provided that 
the outputs are all sampled at a suitably chosen single sampling rate 1/To, and that the inputs are 
updated at the rates Ni/T 0 , . . . , N m /T 0 , where the Nj are certain positive integers. 
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The problem with pole placement for multirate systems is the same as with pole placement for 
single-rate systems: how to determine where the closed-loop poles should be placed? It is a 
particularly difficult problem in the muitirate case because multirate systems are periodically time 
varying (2,81. The periodicity of multirate systems implies that their eigenstructure can only be 
defined based on their (time-invariant) BTP-to-BTP dynamics. Determining desirable closed-loop 
poles for a muitirate system is typically difficult because the BTP of a multirate system is typically 
much longer than the characteristic times of many of its faster dynamics. 

Singular perturbation control law synthesis methods (9,10,11,12,13,14,151 were first developed 
for continuous-time control systems to take advantage of the multiple-time-scale dynamics that 
often occur in control systems. It would seem that an extension to multirate sampied-data systems 
should follow naturally, given that a principal motivation for multirate sampling has always been 
to take advantage of those same multiple time scales, but that has not been the case in practice. 

The problem is the singular perturbation method’s inherent dependence on a coordinate 
transformation to separate the full control law synthesis problem into two (or more) dynamically 
decoupled control law synthesis problems of different time scales. Such a coordinate 
transformation is the first step in control law synthesis by the singular perturbation method. The 
state coordinates are easily decoupled because they represent only the plant s internal dynamics. 

The input and ouput coordinates cannot be so manipulated because they represent the plant's 
external sensor and actuator signals. Consequently, during the second control law synthesis step, 
when the control laws for the different-time-scale state vector components are synthesized 
separately, every control input vector element and every sensor ouput vector element remains 
coupled to every state coordinate so that, just as with successive loop closures, all dynamic coupling 
between control loops cannot be accounted for. 

Various schemes have been developed to circumvent this difficulty. None have been 
completely successful. In Ref. 13, for example, a state feedback control law is synthesized by the 
singular perturbation method, and the lack of a completely decoupling transformation gives rise to 
a requirement for the slow component of the plant state vector to be estimated between slow- 
sampler updates, and a requirement for every control input to be updated at every 
sampling/update instant. 
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The advantage of the LQG method 12,16,17,181 for multirate sampied-data control law synthesis 
is that the control laws for all control loops are synthesized simultaneously, taking into account all 
dynamic coupling between control loops. The disadvantages are the same as with the LQG method 
for continuous-time control law synthesis: that practical performance and stability robustness 
objectives are often difficult to achieve via the minimization of a quadratic performance index, and 
that the resulting control laws are often unnecessarily complex. LQG control laws are even less 
desirable in the multirate as compared to the single-rate case because multirate Kalman filter and 
LQR state feedback gains are periodically time-varying [21. In short, LQG multirate sampied-data 
control laws can provide a useful benchmark for performance comparisions, but they are not 
practical for applications. 

Parameter optimization methods [2,191 for multirate sampied-data control law synthesis 
combine the principal advantages of the LQG and successive loop closures synthesis methods. They 
allow the synthesis of multirate sampied-data control laws of practical structure, and 
simultaneously account for all dynamic coupling between control loops. The typical parameter 
optimization method requires that the control law structure and its parameters to be optimized be 
prescribed. A numerical search is used to determine values for those parameters such that a 
performance index is minimized, possibly subject to constraints on those parameters. The 
disadvantage of parameter optimization methods is that they inevitably require a numerical search 
to determine the control law parameters. 

A new parameter optimization method for synthesizing multirate sampied-data control laws is 
described in Sec. Ill of this paper. It is the second generation of the method described in Refs. 2 and 
20. Unlike its predecessor, which accomodates only partial state feedback control laws, this new 
method accomodates a general, dynamic, multiple-input multiple-output control law structure. 
This new control lawstructure is described in Sec. II of this paper. Section IV describes an 
application of this new method to a design problem involving a two-link robot arm model. 
Conclusions are given in Sec. V. 


II. CONTROL LAW STRUCTURE 
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This section describes the muitirate sampled-data control law structure in Fig. 2. In Fig. 2, y is 
the noise-free, continuous-time sensor signal, v is the discrete-time sensor noise signal, and u is 
the continuous-time control signal. The one sampler in Fig. 2 operates at the sampling rate 1/T, 
where T is the STP of the system s complete sampiing/update schedule. The Delay blocks are one- 
STP delavs. The ZOH block is a zero-order hold. 

The sensor sample-and-hold dynamics are represented by 

v(m,n+l) = (1 - Sy(n)] y(m,n) + Sy(n) y(m,n) 

where v is the sensor signal hold state vector. The matrix Sy(n) is the sensor switching matrix for 
the (n+ l)th STP. We define a switching matrix as a diagonal matrix with 1 or 0 at every diagonal 
position. If the ith diagonal element of S y (n) is 1, the continuous-time signal from the ith sensor is 
sampled at the start of the (n+ l)th STP of every BTP and that sampled value is immediately stored 
as the ith element of y; otherwise, the same element of y is held at those instants. The key point is 
that y always contains the most recent sampled sensor data. 

The processor dynamics are represented by 

z(m,n+l) = (I - S z (n)]z(m,n) + S z (n) (A z (n) z(m,n) 

+ B z (n) {[I - Sy(n)ly(m,n) + Sy(n) y(m,n)}| (2 

u(m,n) = C z (n) z(m,n) 

+ D z (n) ([I - Sy(n)l y(m,n) + Sy(n) y(m,n)j (3 

where z is the processor state vector, and <i is the processor output vector. The matrix ^(n) is the 
processor state switching matrix. If the ith diagonal element of S z (n) is 1, the ith processor state is 
updated at the start of the (n+l)th STP of every BTP; otherwise, the same element of z is held at 
those instants. The matrices A z (n), B z (n), C z (n), and D z (n) are the processor state model matrices, 
whose determination constitutes the control law synthesis problem. Note that a nonzero D z (n) 
results in direct feedthrough of sensor data to u(m,n). 
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The control signal update-and-hold dyamics are represented by 


uirn,n+l) = 1 1 - S, d (n)l u(m,n) + Su(n) u (m,n) 


( 4 ) 


where u is the control signal hold state vector. The matrix Su(n) is the control signal switching 
matrix. If the ith diagonal element of S u (n) is 1, the ith element of u is updated at the start of the 
! n + i)th STP of every BTP; otherwise the same element of u is held at those instants. 


Finally, the continuous-time control signal u is generated by 

u(t) = [1 - S^nll u(m,n) + S u (n) u(m,n) ^ 

for all t on {(mP + n)T, (mP + n + 1)T). 

The advantage of the control law structure of (1) through (5) is that it can be used to represent 
virtually any sampled-data control law structure of practical interest. Its form, however, is not 
standard. Straightforward algebra, applied to (1) through (5), yields the following more standard 

form: 


c(m,n+l) = A c (n) c(m,n) + B c (n) y(m,n) 
u(m,n) = C c (n) c(m,n) + D c (n) y(m,n) 


where 


c(m,n) = (z(m,n) u(m,n) y(m,n)l T 


A c (n) = 


(I - S z {n)l + S z (n) A z (n) 
0 

Su(n) C z (n) 


S z (n) B z (n) (I-Sy(n)l 

I - Sy(n) 

S u (n) D z (n) (I - Sy(n)l 


0 

0 

I - S u (n) 


( 8 ) 


(9) 
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B c (n) = 


( 10 ) 


~S z (n) B z (n) Sy(n) - 
Sy(n) 

_S u (n) D z (n) Sy(n)_ 

C c (n) = (SuCn) C z (n) S u (n) D z (n) [I - S y (n)] I - S^n)] 
D c (n) = [Sy(n) D z (n) Syln)] 


(ID 

( 12 ) 


with 


u(t) = u(m,n) v 

tor all t on ((mP + n)T, (mP + n + 1)T). 

III. PARAMETER OPTIMIZATION METHOD 

This section describes a parameter optimization control law synthesis method for the control 
law structure of Sec. II. It is a generalization of the similar method for state feedback control laws 
described in Refs. 2 and 20, and incorporates also the multiple-plant-condition design for 
robustness ideas of Ref. 21. The approach involves a numerical search to determine the processor 
matrices, A z (n), B z (n), C z (n), and D z z(n), for n=0, . . . , P-1, such that a quadratic performance index 
is minimized. That approach has been criticized in the past because of (1) the difficulties of 
achieving practical performance and stability robustness objectives via the minimization of a 
quadratic performance index, and (2) difficulties related to the convergence of the numerical search. 

The proposed method addresses those criticisms in several ways. First, to enable synthesis for 
robustness to plant parameter variations, the performance index is defined over multiple plant 
conditions. This simple idea has been a key to the success of the popular Sandy [21,22,23,24,25,261 
algorithm for synthesizing robust continuous-time control laws. Second, to improve the 
convergence of the numerical search, the performance index and its gradients with respect to the 
control law parameters are calculated exactly, at every iteration, using closed-form expressions. 
Third, so that a stabilizing initial guess for the control law is not required, and to eliminate 
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problems with destabilizing control laws encountered during the search, a finite-time performance 
index is used. Finally, to lessen the difficulties of achieving practical performance and stability 
robustness objectives via the minimization of a quadratic performance index, linear and nonlinear 
constraints can be imposed on the control law parameters. 


The continuous-time plant dynamics at plant condition i are assumed to be represented by: 


pU) (0 = a" )- 


p (i) (t) + Bi l !u (i) (t) + Bl'!,w a) (t) 


(14) 


F 1 


pw 


v«>. 


(o = e 

r 


p U) (t) 


(15) 


where p a) is the piant state vector, u (i) is the control input vector, y (,) is the sensor output vector, 
and w (l) is a stationary, zero mean, gaussian white noise input vector of known power spectral 
density. 


The performance index is assumed to be 

N, 

J = 


_ N P 

i=l 


Ll? 

'p u> (o' 

T 

Q (i) 0 

p (i) (t)' 

dtr 

I 2 " i 



. 0 R (i) . 


Ul 


(16) 


where Np is the number of plant conditions; E is the expected value operator; tf is the final time 
and is a multiple of the BTP of the system s complete sampling/update schedule; and Q a) and R (i 
are the state and control weighting matrices for the ith plant condition and are non-negative 
definite matrices. 

Based upon the description of the continuous-time plant dynamics in (14) and (15), a complete 
description of the complete system s sampling/ update schedule, the performance index in (16), and 
the control law in (6) through (13), closed-form expressions for the performance index J and for its 
gradients with respect to the processor matrices A z (n), B z (n), Qfn), and D z (n), for n=0, . . . , P-1, are 
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derived in Refs. 19 and 27. Those derivations and the resulting closed-form expressions are 
lengthy, and will not be repeated here. The key points are that the resulting expressions are closed- 
form, and that the number of computations required for their evaluation is independent of tf. The 
single restriction for those expressions to be valid is that the state transition matrix for the BTP-to- 
BTP closed-loop system must be diagonalizable (19). That is not a serious restriction because that 
matrix is rarely nondiagonaiizable in practice. 

Thus far nothing has been said about synthesizing other than periodically time-varying control 
laws. To that end, the performance index and gradient derivations in Refs. 19 and 27 assume that 
the processor matrices are constrained to satisfy 

' D z (n) C,(n) 1 M-l 

I I = Z ct(n,r) 

j_ B z (n) A z (n) J r=0 

with Me (1, . . . , P), and with the a functions constrained to satisfy 

(1 if p = q 
ct(n,p) a(n,q) = j 0 ifp * q 


D z (r) C 7 (r) 

L B z (r) A z (r) J 


(17) 


Equations (17) and (18) constrain the number of different sets of processor matrices to M. The 
function a(r,n) determines which set of processor matrices is active at the (n+ l)th STP. Equation 
(18) guarantees that only one set of processor matrices is active per STP. 

Based on the description of the continuous-time plant dynamics in (14) and (15), a complete 
description of the complete system's sampling/ update schedule, the performance index in (16), the 
control law in (6) through (13), the constraint relations in (17) and (18), and the closed-form 
expressions for the performance index J, and for its gradients with respect to the processor matrices 
A z (r), B z (r), C z (r), and D z (r), for r=0, . . . , M-l, in Refs. 19 and 27, we have developed a computer 
algorithm to numerically determine a set of processor matrices that minimizes J. A numerical 
search is used to determine the processor matrices A z (r), B z (r), C z (r), and D z (r), for r = 0, . . . , M-l, 
given an initial guess for those matrices. The NPSOL nonlinear programming algorithm is used 
for the numerical search. NPSOL [281 is a powerful nonlinear programming package with good 
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convergence properties as a result of its use of exact performance index and gradient evaluations at 
every iteration. In addition, NPSOL accomodates linear and/or nonlinear constraints on the 
independent variables. This means that linear and/or nonlinear constraints on the control law 
parameters can be combined with the usual performance index minimization objectives to achieve 
practical performance and stability robustness objectives. 


Additional important features of this new synthesis algorithm include automatic discretization 
of the continuous-time plant model and of the continuous-time performance index (27). These are 
important design features because they effectively decouple the sampling/ update rates selection 
problem from the problem of determining a suitable performance index. This means that the 
performance index can be determined first, based on a continuous-time design, and that this new 
algorithm can then be used to determine a multirate sampled-data design that minimizes the same 
performance index. 


In summary, the inputs required to apply this new synthesis algorithm are the following: 

• A state model description of the continuous-time plant dynamics at each of the Np 

plant conditions. 

• State and control weighting matrices for the performance index at each of the Np plant 

conditions. 

• The final time tf for the performance index. 

• The power spectral density of the continuous-time white process noise at each of the 

Np plant conditions. 

• A complete description of the complete system's sampling/update schedule. 

• The integer M and the a functions that constrain the periodicity of the processor 

matrices via (17) and (18). 
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• The desired dynamic order and structure for the processor matrices. 

• The covariance matrix for the discrete-time sensor noise at each of the N p plant 

conditions. 

• A complete description of all linear and/or nonlinear constraints to be imposed on the 

elements of the processor matrices. 

• An initial guess for the processor matrices. 

A disadvantage of most parameter optimization control law synthesis methods is that they 
require a stabilizing initial guess tor the control law. That is not the case with this method because 
of its finite-time performance index. The finite time ensures that the performance index and its 
gradients will be finite whether or not the closed-loop system is stable. A disadvantage of the finite 
time is that a steady-state solution, i.e., for 

T A . , (19) 

Jss= llm J 

tf *4 


cannot be obtained directly. A steady-state solution is easily obtained, in practice, however, by 
choosing a finite time tf that is large compared to the characteristic times of all of the closed-loop 
system's poles. Because the number of computations required to evaluate the performance index 
and gradient expressions of Refs. 19 and 27 does not depend upon t f , this can be done without 
penalty in terms of the computation time for the numerical search. 

In practice, because digital computers cannot store arbitrarily large finite numbers, a steady-state 
solution usually cannot be obtained by simply initially setting tf to a large value. Instead, it is 
usually necessary to complete first (i.e., when the current best guess for the control law parameters 
is poor) an optimization for a small tf, and to then re-optimize, for larger and larger tf, until tf gets 
to be large compared to the characteristic times of all closed-loop poles. 


41 


The final key issue regarding this new synthesis algorithm concerns the requirement that the 
structure of the processor matrices be specified. The key point is that the imposed structure should 
guarantee that the free parameters for the numerical search constitute an independent set with 
respect to the control law s input-output dynamics. When the processor dynamics in (2) and (3) are 
considered, with the constraints in (17) and (18) in effect it is straightforward to see that the 
complete set of the elements of A z (r), B z (r), C z (r) and D z (r), for r = 0, . . . , M-l, do not constitute 
such an independent set because, for example, an arbitrary change in one element of B z (0) can be 
compensated for by changes to the elemets of C z (0), and to the other elements of B z (0), such that the 
processor s input-output dynamics are unchanged. 

Thus, additional structure, or, equivalently, additional constraints, must be imposed on the 
elements of the processor matrices to guarantee that the free parameters for the numerical search 
constitute an independent set with respect to the control law's input-output dynamics. In practice, 
a suitable set of such constraints can frequently be determined based on "classical" control law 
structures (e.g., combinations of lead and lag compensators and notch filters). 


More generally structured control law can, of course, also be accomodated. What constitutes an 
optimal structure for the processor matrices for the general case is a topic of current research. We 
have successfully applied the following structure (shown for the n-is-even case) for the particular 
case where the constraints in (17) and (18) are applied with M = 1 (the time-invariant case): 


A z (0) = block diag 


0 

2 2 
-Oj-tOj 


1 

- 20 , 



( 20 ) 


B z (0) = 


C z (0) = 


b,i ^Tm 

b n i bim _ 

1 • • • 1 

C 21 ' ‘ ‘ c 2n 

c pt c f*n — 


( 21 ) 


(22) 
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D z (0) = 




( 23 ) 


L dp] dpm J 

The appendix shows that the (m+ p)n+ mp ctj, co ( , bjj, Cjj and d S j parameters of this structure 
constitute an independent set with respect to the control law's input-output dynamics (for any 
non-tnvial sampling/update schedule) provided that no eigenvalues Oj ± jco; of the processor 

dynamics are repeated. 


IV. EXAMPLES 

This section describes the design of a tip position control system for a planar two-link robot 
arm. The robot arm system is shown in Fig. 3. The first link is long and massive, for large-scale 
slewing motions. The second is relatively short and lightweight, so that high-bandwidth control of 
the arm's tip position can be achieved using a relatively small motor at the second joint. The pin 
joint, rotational spring and rotational damper at the midpoint of the first link models flexibility in 
that link. The second link is assumed to be rigid. The motor torques Tj and T 2 are the control 
inputs, and it is assumed that only the joint angle 0 and the tip position 5 are measured. The 
linearized dynamical equations for this system for small e - 0 and small (j) - e are easily derived. 
The spring constant (k) and damping coefficient (b) values (in Fig. 3) were chosen based on that 
model to achieve 1 percent damping and a 10 Hz natural frequency for the open-loop vibration 

mode. 


Figures 4 through 6 show the closed-loop arm responses, based on the linearized arm 
dynamics, to a step change in the commanded tip position with nine different control laws. The tip 
position (8) responses are shown in Fig. 4. The simultaneous control torque CT| and T 2 ) responses 
are shown in Figs. 5 and 6. The nine different control laws are briefly described as follows: 


LQR Analog: The continuous-time LQR (full-state-feedback) control law that minimizes 


J = 


t H ® 2tf 


p 2 [02 + (p6) 2 ] + 


2 i _\2 

lmax 



(24) 
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with 


(3 = 32 rad - ' 1 


(25) 


p = 14 rad/m 


(26) 


T 2m ax = 0.00335 N m 


(27) 


T lm ax/T 2m ax = 8 UO/ 

The T-> max value is the To torque that achieves® = 2jc rad/sec in 1 sec with 0(t) s e(t) s 0 andd(O) = 0. 
The T 1max /T 2max value represents a typical ratio of peak motor torques at the respective joints. The 
p and p values were chosen by trial and error to achieve the closed-loop poles in Table 1. Note that 
the ratio of the characteristic frequencies of the rigid-body closed-loop pole pairs is eight; and that 
the characteristic frequency of the faster closed-loop rigid-body pole pair is a factor of five less than 
the characteristic frequency of the closed-loop vibration mode. 


Third-Order Analog Successive Loop Closures: The third-order, continuous-time, successive loop 
closures control law in Fig. 7, which consists of a single lead compensator in the 0-to-T! loop and 
twin, cascaded lead compensators in the 8-to-T 2 loop. The closed-loop poles for this design are in 
Table 2. Note that the rigid body and vibration mode closed-loop poles match those of the LQR 
Analog design. 

Third-Order Multirate Tustin: A muitirate sampled-data approximation to the Third-Order 
Analog Successive Loop Closures design obtained via Tustin's approximations of the continuous- 
time transfer functions in Fig. 7. The sampling/update rates (in sampies/updates per second) of the 
0-to-T! and 8-to-T 2 loops are eight times the characteristic frequencies (in cycles per second) of the 
slow and fast, respectively, rigid body closed-loop pole pairs from the Third-Order Analog 
Successive Loop Closures design. 
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Optimized Third-Order Multirate Tustin: The same as the Third-Order Multirate Tustin design, 
but with the lead compensator gain, zero, and pole locations optimized, by the paramether 
optimization control law synthesis method of Sec. Ill, to minimize the same performance index as 
in the LQR Analog design. To synthesize this control law, continuous-time process noise and 
discrete-time sensor noise inputs were added to the robot arm model. The former were taken to be 
white noise disturbance torques and w 2 , coincident with the respective control torques, with 


f wi(t)l I 4.9x10 

El wLri h w l (t) w 2 (t) l} = n 


0 


8(t — t) 


(29) 


where 6 is the Dirac delta function. The latter were taken to be stationary, purely random 
sequences. v 1 and vi, for the 9 and 5 measurements, respectively, with 

E ([v2(m,n)J [v l (m ' n) V 2 (m ' n)1) = 

Multirate Third-Order; The same as the Optimized Third-Order Multirate Tustin design, but using 
the third-order, generalized, time-invariant structure in (20) through (23) for the processor 
matrices. Just as with the Optimized Third-Order Multirate Tustin design, two of the processor 
states are updated at the faster sampling/ update rate, and the third is updated at the slower 
sampling/ update rate. 

Multirate Second-Order: The same as the Multirate Third-Order design, but using the second- 
order, generalized, time-invariant structure in (20) through (23) for the processor matrices. One of 
the processor states is updated at the faster sampling/update rate, and the other is updated at the 
slower sampling/update rate. 

Multirate First-Order: The same as the Multirate Third-Order design, but using the first-order, 
generalized, time-invariant structure of (20) through (23) for the processor matrices. The one 
processor state is updated at the faster sampling/ update rate. 


8.1 xlO -3 

0 


1x10" 


(30) 
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Single-Rate Third Order: The same as the Multirate Third-Order design, but single-rate, with the 
single sampling/update rate chosen to yield the same number of real-time computations per unit 
time as the Multirate Third-Order design. 

Analog Third-Order: The continuous-time equivalent to the Multirate and Single-Rate Third- 
Order designs. The processor matrices have the same structure as in the Multirate and Single-Rate 
Third-Order designs. The control law was synthesized using the Sandy algorithm (211 to minimize 
the same performance index as in the Multirate and Single-Rate Third-Order designs. 

The LQR Analog responses in Figs. 4a, 5a and 6a constitute the optimal responses for the 
performance index in (24), assuming full state feedback, no process or sensor noise, and infinitely 
fast sampling. The Third-Order Analog Successive Loop Closures responses in the same figures 
have low tip position overshoot, but include also a relatively large contribution from the vibration 
mode (see especially Figs. 5a and 6a). 

The Third-Order Multirate Tustin responses in Figs. 4a, 5a and 6 are unacceptable. This is 
somewhat suprising, but not totally unexpected given the (low) factor-of-eight sampling/update 
rate-to-characteristic frequency ratio for this design. 

The Optimized Third-Order Multirate Tustin responses in Figs. 4a, 5a and 6a are acceptable, and 
demonstrate that the parameter optimization control law synthesis algorithm of Sec. Ill can be used 
to optimize the parameters of classically-structured control laws. 


The Multirate Third-Order, Second-Order and First-Order responses in Figs. 4b, 5b and 6b 
demonstrate that the same parameter optimization control law synthesis algorithm can be used to 
synthesize multirate sampled-data control laws having a prescribed dynamic order and a 
prescribed, but general, structure, with apriori specified sampling/update rates for all sensors, 
processor states, and control inputs. 

The Single-Rate Third-Order and Analog Third-Order responses in the same figures put the 
multirate responses in perspective. The Single-Rate Third-Order control law is the single-rate 
equivalent to the Multirate Third-Order control law because it (1) was synthesized to minimize the 
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same penormance index, using the same process and sensor noise characteristics: and <» requires 
the same number of computations per unit time for real-time operation. 

The Analog Third-Order responses in the same figures are the responses that would have been 
obtained with the Multirate Third-Order control law and Singe-Rate Third-Order control laws if 
sampling and update rates were not an issue, and very fast sampling and update rates were 

everywhere used. 


V. CONCLUSIONS 

With the possible exception of successive loop closures, the multirate sampled-data control law 
synthesis methods available today fail to provide the designer with sufficient flexibility to prescribe 
sensor sampling rates and processor state and control input update rates. A new parameter- 
optimization-based method for synthesizing multirate sampled-data control laws of arbitrary 
dynamic order that provides that flexibility is described in this paper. This new method, described 
in Sec. in, determines, by numerical optimization, the free parameters of the general purpose 
multiple-input, multiple-output, sampled-data control law structure in Fig. 2, to minimize a 
quadratic performance index, possibly subject to linear and/or nonlinear constraints on those 
parameters. A stabilizing initial guess for the control law is not required because the performance 
index is finite-time. To enable the synthesis of robust control laws, the performance index can be 
defined over multiple plant conditions. 

An application of this new method to the design of a tip position control system for a sixth- 
order, two-link robot arm was described. Multirate sampled-data control laws of various dynamic 
orders synthesized by various methods were compared to confirm that the new synthesis method 
can be used to synthesize multirate sampled-data control laws having a prescribed dynamic order 
and structure, with apriori specified sampling/ update rates for all sensors, processor states, and 

control inputs. 
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APPENDIX A. SUPPORTING MATHEMATICS 

Consider the control law in (1) through (5). Suppose that the constraints in (17) and (18), with 
\1=1, are in effect, so that the processor matrices are constrained to be time-invariant. Suppose that 

the processor matrices. A z (0), B z (0), C z (0) and D z (0), are further constrained to have the forms in 
20) through (23). Finally, to guarantee a nontrivial sampling/ update schedule, suppose that the 
sensor, processor state, and actuator switching matrices satisfy 


r p-i 1 

det| £s y (n) j * o 

L n=0 

r p-i i 

detl S 2 (n) I * 0 
L n=0 J 

r p-i 1 

detl £ S u (n) j * 0 
n=0 


(31) 


(32) 


(33) 


We will show that the (m+p)n+pm ctj, co , bjj, Cjj and d t j elements of the control law then constitute 

an independent set with respect to that control law's input-output dynamics if and only if A z (0) has 
no repeated eigenvalues. 

We begin by noting that, with (31), (32) and (33) in effect, it is straightforward to see that the 
independence in question does not depend whatsoever on the sensor, processor state, or actuator 
switching matrices. Therefore we consider only the special case where S y (n), S z (n) and S u (n), for 
n=0, ... P-1, are identity matrices. The control law then reduces to 


z(m,n+l) = A z (0) z(m,n) + B z (0) y(m,n) 


(34) 


u(m,n) = C z (0) z(m,n) + D z (0) y(m,n) 
where u(m,n) is defined in (13). 

Consider first the control law 
z(k+l) = A z(k) + By(k) 
u(k) = C z(k) + D y(k) 

where 



'zf 


_ yi" 


-u,l 

z = 

-Zn- 

li 

-Ym- 

u = 

1 

"O 

1 


A = 


0 

0 


0 0 
0 


o K J 


(35) 


(36) 

(37) 


(38) 


(39) 
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and B, C and D have the forms in (21) through (23). So that the control law s impulse response will 
be purely real, the and the corresponding columns of C and rows of B must be either real, or 
must occur in complex-conjugate pairs, and D must be pureiv real. 

Lemma: The (m+p)n+pm k\, b^, c^, and d jj parameters (counting a real element as one parameter, 
and a complex conjugate pair of elements as two parameters) of the control law in (36) and (37) 
constitute an independent set with respect to that control law s input-output dynamics if and only 


if X,*Xj, for i*j. 


Proof: Consider the related control law 


z(k+l) = A z(k) + B v(k) 

(40) 

u(k) = C z(k) 

(41) 


with 




Its inputnautput dynamics are represented by 


n 

H(z)=C(zI-Ar 1 B = 

i=l 



QB t 

z-Xj 


(42) 


(43) 


where Q is the i th column of C, and Bj is the ith row of B. The bjj and Cjj parameters of this 
control law are dependent with respect the control law's input-output dynamics if and only if, for 

an arbitrary change in one, the others can be changed so that H(z) is unchanged. 

Case 1: No repeated X; 's. 

From (41), for the case of no repeated Xj 's, it is straightforward to see that when one of the Xj's is 
changed by an arbitrary amount, it will not be possible to change the remaining X,, bjj and c ( j 
elements so that H(z) is unchanged. But when bjj is multiplied by a nonzero but otherwise 
arbitrary a, we can multiply the remaining elements of B } by a, and divide Cj by a, so that H(z) is 

unchanged. Thus, the X,, bjj and Cjj parameters of the control law in (40) and (41) are dependent 
with respect to that control law s input-output dynamics. 

If, however, one element of every column of C is fixed, as is the case in the C matrix of (22), it is 
similarly straightforward to see that it will not be possible to compensate for an arbitrary change in 

any X;, bjj or Cjj element by changing the remaining Xj, bjj and Cjj elements so that H(z) is 

unchanged. Thus, for the case of no repeated Xj's, with one element of every column of C fixed. 
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the A.,, bjj and parameters of the control law in (40) and (41) constitute an independent set with 
respect to that control law s input-output dynamics. 


Case 2: Repeated A.,'s. 

Consider again the control law in (40) and (41). but this time suppose that A.i=A 2 . That control law's 
input-output dynamics are represented by 


. QBt+CoBo 
H(z) = C (z I-A) B = : 

Z-A-i 



C.Bj 

z-A., 


(44) 


From 142), it is straightforward to see that, with or without one element of every column of C fixed, 
the remaining bjj and c,j elements of C v B,, C 2 , and B 2 can be changed to compensate for an 
arbitrary change in any one element of C-j, B,,C 2 , or B^ so that H(z) is unchanged. Thus, for the 

case ot repeated AT s, with or without one element of every column of C fixed, the A.;, bjj and c,j 
parameters ot the control law of (40) and (41) are dependent with respect to that control law’s input- 
output dynamics. 

General Case: We conclude that the (m+p)n A* bjj, and c,j parameters of the control law in (36) and 
(37), with D= 0, constitute an independent set with respect to that control laws input-output 
dynamics if and only if A.,*A.j, for i*j. A nonzero D matrix simply adds pm parameters to that set. 

Theorem: For the control law in (1) through (5); with the constraints in (17) and (18), with M=l, in 

effect, so that the processor matrices are constrained to be time-invariant; with A z (0), B z (0), C^O) 

and D z (0) further constrained to have the forms in (20) through (23); and assuming that the sensor, 

processor state, and actuator switching matrices satisfy (31) through (33); the (m+p)n+pm o,, (0,, bjj, 
c^ and djj elements of that control law constitute an independent set with respect to that control 

law's input-output dynamics if and only if A z (0) has no repeated eigenvalues. 

Proof: The control law in (34) and (35) has the same number of free parameters as the control law 
in (36) and (37), and the two are related by the similarity transformation x= M z , where 


M = block diag 


1 I" 


2j 0), 


L -<T,+jO), 


1 

-1 


1, . . . , n/2) 


with 0 , = RefA.,) and to, = ImKA.,). 


(45) 
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Closed-Loop 

Pole 

Damping 

Ratio 

Characteristic 

Frequency 

Rigid Body 

-1.10 ± 1.10 

0.71 

0.25 Hz 

Rigid Body 

-8.81 ±; 8.83 

0.71 

2.0 Hz 

Vibration Mode 

-0.649 ±j 62.8 

0.01 

10 Hz 


Table 1 LQR Analog Design Closed-Loop Poles 



Closed-Loop 

Pole 

Damping 

Ratio 

Characteristic 

Frequency 

Rigid Body 

-1. 10 ±j 1.11 

0.71 

0.25 Hz 

Rigid Body 

-8.88 ±; 8.84 

0.71 

2.0 Hz 

Vibration Mode 

-1.35 ±j 63.9 

0.02 

10 Hz 

Compensator 

-10.5 

- 

1.7 Hz 

Compensator 

-33.2 ±; 34.0 

0.70 

7.6 Hz 


Table 2 Third-Order Analog Successive Loop Closures Design Closed-Loop Poles 
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Time Lines for Sampling/Update Activities: 



0 4T 8T 12T 16T 20T 24T 



Fig. 1 Example Multirate Sampling/Update Schedule 
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Fig. 2 Multirate Sampled-Data Control Law Structure 
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Parameters: Mass Length 

L\ 0.5 kg 0.5 m 

Lo 0.5 kg 0.5 m k = 37.33 N/rad 

L 3 0.04 kg 0.2 m b = 0.012 N- s/m 

The natural frequency of the vibration mode is 10 hz. 

Inputs: Torques T i andT 2 

Outputs: 0 and 8 


Fig. 3 Two-Link Robot Arm 
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LQR Analog 

Third-Order Analog Successive Loop Closures 
Third-Order Multirate Tustin 
Optimized Third-Order Multirate Tustin 


Seconds 


Fig. 4a Tip Position (8) Responses to Tip Position Step Command 


Analog Third-Order 
Multirate First-Order 
Multirate Second-Order 
Multirate Third-Order 
Single-Rate Third-Order 



Seconds 


Fig. 4b Tip Position (8) Responses to Tip Position Step Command. 
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Fig. 5a Control Torque T\ Responses to Tip Position Step Command 
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Fig. 5b Control Torque T\ Responses to Up PosmoaStep Command- 
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i 1.5 2 2.5 3 

Fig. 6a Control Torque T 2 Responses to Tip Position Step Command 


3.5 4 

Seconds 



Fig. 6b Control Torque T 2 Responses to Tip Position Step Command 
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Introduction 

In many cases a muitirate compensator can provide better per f or m ance than a single rate 
compensator requiring the same number of computations. Berg, for example, was able to 
reduce the steady state RMS response of states and controls for a simple mass-spnne-mass 
system neariv 20% by using a muitirate compensator over a singie rate compensator. 
Numerous otner examples have been provided in the literature by Bere [3]-[5], Amitt (1]- 
[2], ana Yang [ 17] . While multirate compensators can provide improved performance over 
singie rate compensators, they are also, in general, more complicated to desien. 

The complexity of multirate compensators stems from the fact tnat they are by nature 
:ime varying, periodically time varying for any practical application. Not only must 
designers cnoose muitmie sampiing/update rates for the compensator, cut mev must also 
determine the parameter values tor a time varying compensator. 

We wiil consider three methods which can be used to design muitirate compensators: 
Multirate LQG [ 1 ]; Generalized Algorithm for Multirate Synthesis iGAMS) [ 17]; and our 
new method presented in this paper. All three of these methods synthesize linear 
compensators by minimizing a cost function quadratic in the states and controls. 

Multirate LQG is the multirate equivalent of singie rate LQG. The muitirate LQG 
problem is straightforward to soive because the equations governing the solution are similar 
to those for the singie rate case so that most methods used to soive the sineie rate LQG 
problem can be applied to muitirate LQG. The disadvantage of muitirate LQG is that it 
results in a lull order compensator which has periodically time varvine gains. 

GAMS was developed by Yang to overcome many of the short comings of multirate 
LQG. Yang s algorithm can synthesize reduced order multirate compensators with or 
without time varying gains by using a numerical gradient type search to find the optimum 
compensator values. His algorithm uses a finite time cost function in its problem 
formulation, unlike muitirate and single rate LQG which use an infinite time cost function. 
By using a finite time cost function Yang’s algorithm eliminates the numerical problem that 
arises when a destabilizing compensator is encountered during the numerical search since 
the value of the cost function is infinite at infinite time when the closed loop system is 
unstable. The primary benefit is that the designer does not have to provide the numerical 
search algorithm with a initial stabilizing guess. Unfortunately, using a finite time cost 

function greatly complicates the equations which govern the soiudon, making them 
laborious to solve. 
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We have been using Yang's algorithm to soive muitirate problems for the iast two years 
and have found that the flexibility afforded by a finite time cost function is not needed to 
soive most muitirate oroolems. For this reason we have developed a third method for 
synthesizing muitirate compensators. 

Our new algorithm is similar to Yang's (GAMS) in that it can be used to synthesize 
reduced order muitirate compensators with or without time varying gams: but uniike 
Yang s aiaonthm we use an infinite time cost function in the problem formulation. The 
advantage of this new approach is that it resuits in a set of relatively simple governing 
equauons which are closeiv related to those for the single rate case. 

In this paper we will present our new algorithm. In Section I we discuss the structure 
of the general muitirate compensator. In Section H we develop the equauons which govern 
the soiudon of the reduced order muitirate compensator. Section III contains a brief 
discussion of how we implemented our algorithm. In Section IV we present a simple 
muitirate example followed by some concluding remarks in Secuon V. 

L The General Multirate Compensator 

Before deriving the equations governing a reduced order muitirate compensator, we 
will first present the structure for a general muitirate compensator. We restrict our 
discussion for now to compensators with time invariant gains and sampiing/update rates 
whose ranos are rational numbers. 

A general muitirate compensator is shown in Figure 1. Each input (y), output (u), and 
state tz) is sampied/updated at a rate which, in general, represents the desired bandwidth of 
the input or output with which it is associated, y is the value of y currently available to the 
digital processor from the zero order hold: while u is the current output from the digital 
processor which is held with a zero order hold to form the output u. When the 
sampiing/update rates have ratios which are rational numbers the sampling schedule is 
periodically time varying. The greatest common divisor of all the sampling/update periods 
is the shortest time period (STP); the least common multiple of all the sampling/update 
periods is the basic time period (BTP) (see Figure 2). 

The state equations for the multirate compensator pictured in Figure 1 are: 

| - \ [l*s z ,icj+s z , k A s z<k B[l-s y ,it] 

y j = ' 0 [J-Sy.kj 

' U ^ L S u*C S u .kD[l-Sy.k] 
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. _ _ . .A\\ - - - 

u k- ; _Su,kC S u .k.Oj_I-Sy,kJ L I Su.kjj y j *'_ s u,kOs v .k| yk (2) 

1 u Ik 

u is a hold state used to model the sampler and zero order hold between ii and u. s v k , s z 
and s u k are switching matnees tor v, z, and u respectively that model the system's 
^ampimg/update activity at the start the k 1 * 1 STP. s* ^ has the form: 

r n o o ••• o 1 

1 i 

0 r 2 0 ••• 0 | 

s.k = : : | 

0 0 0 

0 • ■ • 0 0 j 

j 1 if the j :h (z, y, or u) is sampied/updated 
where ^ = < at the start of the k th STP 

1.0 otherwise 

m^ = the number of states (z) 
nx, = the number of inputs (v) 
nty = the number of outputs tu) 

A more complete discussion of this compensator structure can be found in [ 171. 

Equanons (1) and (2) can be written more compactly as: 

z k+l = ^k z k + ®kyk ^ 

u k = C k z k + D kVk W 

( z l 

where Zk = ; y j 
1 u |k 

Equations (3) and (4) form a single rate periodically time varying system with a 
sampling rate of one STP and a period of one BTF. If N=BTF/STP, then A k = A fc+N , 

B k = B k+N* C k = C k+N» 111(1 D k = D k+N‘ 
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Even though A k , B k , C k and D k are periodically time varying, the multirate 
compensator has ume invariant gams. The periodicity of the muitirate compensator is due 
:o muitirate sampiing/updadng not the compensator gams. In the remainder of this section 
we wiil demonstrate how tne ume invariant compensator gains. A. B. C, and D can be 
separated from the periodic compensator matrices A k , B k , C k and D V: . 

Defme the composite compensator matrix as 


Pk 


i D * 
L B k 


Ck i 

A k j 


(5) 


and factor P k as follows: 
where 


P k ” S lk PS 2k ' r S 2k 


P = 


D C 1 

Lb aj 


Sik = 


Su.k 

0 

0 

Su.k 


0 

s z.k 

0 

0 


c _[* s y,k 0 I“ s v,k 0 I 

s *-| o i d o i 


(6) 

(7) 


( 8 ) 


(9) 


| 0 0 0 I-Su, k , 

s J 0 I-s z . k 0 0 I 

3k Sy,k 0 I-Sv.k 0 I 

L o oo i-s^j 


( 10 ) 


Equation (6) is a key result. It allows us to factor the time invariant compensator gains* 
the unknown parameters we will solve for in Section II, out of the time varying 
compensator. 

It is important to note the difference between P k and P in (6). P k (with a subscript) is a 
periodically rime varying matrix defined by (5). It includes all the information about the 
compensator gains and the sampling/update schedule. P (without a subscript) is a constant 
matrix which contains only the gains for the compensator. P k can be written in terms of P 
and S lk , So k , and S3 k , equation (6), where S lk , S 2k , and S 2k are periodically rime varying 
matrices which contain a description of the sampling scheme. 
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II. DERIVATION OF REDUCED ORDER MULTIRATE COMPENSATOR 

In this section we wili use me results of Secuon i to derive the equauons governing the 
reduced order muidrate compensator. The muitirate problem to be soived is as follows: 
Given: tne system 

Xk+i = Fx k + Gu k + Wwi( ( 1 1) 

y k = Hx k + v k (12) 


where F. G, W and H are obtained by discretizing the analog plant matrices at one STP; 
and where w k and v k are discrete-rime gaussian white noise inputs. 

Find: the muitirate control law with a prescribed order and sampling scnedule, of the 
form ot 1 1 M2), which minimizes tne Quadratic cost funcuon of the form: 


J = 


N 

lim X E i 


k=l 


J|Xk|T| 
|Wlc) i 


Qi 

m t 


M 

02 



(13) 


E is the expected value operator, and the summation from 1 to N accounts for the fact that 
the closed loop system is periodically time varying. A prescribed sampling schedule implies 
that the values of s 2 j., s y k , and s^ are known. Sample schedule selection is discussed in 

[4] and is the subject of future research. 

When we write 1 1)-(2) as (3)-(4) it is easy to see that this problem is essentially a time 
varying feedback problem - a time invariant plant with a periodically time varying 
compensator. What makes this problem difficult is that the time varying compensator has 
an explicit form, that of (l)-(2), in which only certain parameters. A, B, C, and D. can be 
adjusted to mi nimize i. 

To solve the muitirate control problem we cast it into output feedback form and follow a 
derivation similar to Mukhopadhyay ’s for the single rate case (121, [141. Using (l)-(2) 
written as (3)-(4), and (1 1)-(12) we write the output feedback equations: 


I *♦' 1 
1 I 


F 0 

ML 

" ys, 

G 0 

1 “k 1 

u 

W 0 

? 1 

o oj 

Ml 

0 0 

\ Zk + l 1 

1 

0 0 

l v k 1 


II 

- 

H 0 

l z W 

_ 0 0 




I 

Zkj 

. 0 

0. 


wic ( 
Vic I 


(14) 


( 15 ) 
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" k U D» 

Zk+1 | B * 


Ck ',m 

A k j|z k | 


Eauanons ( 14)-( 16) can be wnnen more compactiv as 


x k + i = Fx k + Gu k + WT lk 
v k = Hx k + Vr( k 

u k = p !<yk 


(16) 


(17) 

(18) 
(19) 


It is important to keep in mind that P k in equation (19) corresponds to the P k in equation 
' 5), a periodically time varying matrix which contains ail the informauon about the muitirate 
compensator gams and samoling/update rates. 

The closed loop system is 


Xk+1 — E ck X k t G c kT]k 

(20) 

Fck = F + GP k H 

(21) 

G ck = W + GP k H 

(22) 


The state covariance propagation for this system obeys the following equation: 

Xk + i = F ck X k Fj k + G ck RGc k (23) 

where 

X k = E(x k x{) 

R = Elriklk) 

Equauons t20)-(22) represent a periodically time varying system with a period of N. 
We can generate a single rate system by repeated application of equation (20) over one BTP 
[11]. The single rate system can be written as 

*k+N = F bk*k + GblcHbk (24) 

where 

Fbk = Fc(k+N-l)Fc0c+N-2}Fc(k+N-3V • ’Fck (25) 

Gbk = S.Fc(k+N-nFc(k+N-2V • -FcOc+oGck j 

Fc{k+N-l)Fc(k+N-2)-*-Ec(k+-2)Gc(k+l) | ••• | G c (k+N-1>] 


68 


1991 ACC 


Reduced Order Multirate Comocnsator Synthesis 


7 


Tbk 


% 

Tk+1 


Tk+N-1 


This single rate system has exactly the same values tor x as the periodically time 
varying closed loop system at each BTP. However, the values of x at the intermediate 
STP’s have been lost because x is incremented by N in (241 but only by l in (20). There 
are actually N such single rate systems associated with (20). They can be written as: 

x k+N+i = Fb(k+n x k+i + G*k-H)Tlk + i fori=l,2...,N (27) 

If F bk is stable, then the periodically rime varying system (20) is stable f9]. We can 
.aicuiate the steady state covanance for x using the following Lyapunov eauauons: 

X k = F bk X k Fj k + G bk R b G bk for k=l .2 N (28) 


Rb = 


R 0 • 
0 R- 


0 

0 


0 0 0 R 


Note that the covanance of x is periodic, that is the covanance vanes within one BTP, 
but from BTP to BTP X, ; . = X«. +N . Once we have calculated the steady state covariance for 
x k at any k using (28), we can use (23) to propagate the covanance over the BTP. This 
eliminates the need to solve equation (28) N times. 

Now, using (23) and (13) and the properties of the Trace (Tr) operator we can write the 
cost function for the stabilized system as 

N 

1 = X Tr ([Qi + MP *H + (MP k H) T + (P k H) T Q 2 P k H|X k + (PkV) T Q 2 PkVR) (29) 

k=l 


Adjoin the constraints to the cost J using Lagrange multipliers, to obtain: 

N 

J = X Tr (fQi + !^PkH +(MP k H) T +(P k H) T Q2p k HjX k +(P k V) T Q2PkVR 
+ ^k+i F ck X k Fg k + G ck RGj k - X k JJ 

with X| =X N+ j. 


( 30 ) 
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The necessary conditions tor nunimum J are 


5J .=0. 51 


ax k 


dA k +i 


n • cJ .. 

= 0. ana — = 0 
a P 


(31) 


d2j 


In addition. — - must be positive definite tor a minimum J. 

5P2 

Substituting (30) into (31) and replacing P k with P k =S lk PSa k +S-i k from (6) we obtain: 


dJ 

ax v 


- 0 = Qi + MPvH + (MP k H) T + (P k H) T Q 2 P k H + Fj k A k+1 F ck - A k 
for k=l,2,.„,N with A k = A k+N 


(32) 


oJ 


= 0 = F ck X k Fc k T G ck RGI k - X k+i 


dA k+ i 

for k = 1 ,2 N with X k = X k+N 


(33) 


^=0=2^j k ([(^+G T A k+1 Gjp k [HX k H T +VRV T ] 

+[M T +G T A k+ iF]x k H T js 


T 

2k 


(34) 


Equanons ( 3 2)- (34) are a set of coupled equations, two sets of Lyapunov equations and 
one Riccau eauaaon. which make up the necessary conditions for which P. the multirate 
compensator gain matrices A, B. C, and D in equation il)-(2), wiil minimize the cost 
function J. The values of A, B, C, and D found by solving (,32)-(34) can be substituted 
into ( l)-(2), along with the definition of the sampling schedule. s zJt , s uJc , and s yJc , to form 

the complete time varying multirate compensator. 

To ensure that the compensator gains satisfying (32)-(34) minimize J, we must also 
check that the hessian of J with respect to P is positive definite. In our present algo rithm, 
we do not calculate the hessian explicitly, but use an approximate value calculated by the 
numerical search algorithm discussed in the next section. 

Equations ( 32)-(34) were derived assuming time invariant compensator gains. We can 
easily derive the equations for periodically time varying gams, so that A, 13, C, and D in 
equation ( 1) are periodic. Let 
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with the restriction that A;*n = A;, B,+n = B;, C,*n = C,, and D^n = D. 

'o that 

-ID, Cj ] 

p = Pj = i _ _ ; 06 ) 

L B i A i J 

Then reDiace P with P; and differentiate (30) with respect to Pj to obtain 


^L=0= sy k ![Q 2 +G T A k+l G]PkfHX k H T +VRV T l 

oPj 1 J L ' (37) 

■h M T +G T A k+1 F |X k H T j Sl k for j = i .2....N 

Tnus for every new set of compensator gains we add one new equation of the form of (37). 

Eauauons (32M34) are very similar to the single rate equations. In fact, if we set N=1 
and S , k and equal to the identity matrix we obtain the exact results derived by derived 

by Mukhopadhvay for a single rate case [12]. 

m. Implementation 

In order to determine the reduced order multirate compensator which minimizes the cost 
function J. we need to solve (32)-(34). A numerical gradient-type search algorithm was 
implemented in MATLAB [16] to solve these equations. A flow chan of the algorithm is 
shown in Figure 3. The equations necessary to solve for the Lagrange multipliers, (A.3)— 
(A. 4), are located in the Appendix. 

The algorithm automatically discretizes the analog plant, weighting matrices and 
process noise covariance matrix. Refer to [3] for a discussion of the relevant discretization 
procedures. To ensure that the solution represents a minimum J, the algorithm checks that 
the hessian of J with respect to the free parameters in P is positive definite at the solution 
point. 

Because t32)-(34) are not valid when the closed loop system is unstable, the numerical 
algorithm 1 ) must be provided with an initial stabilizing compensator, and 2) must result in 
a stabilizing compensator at every iteration of the optimization. In our experience, finding 
an initial stabilizing compensator is generally not a problem using successive loop closures 
[4]. To avoid the problem of destabilizing compensators during the iteration process we 
included a check in the algorithm which systematically reduces the step size to ensure the 
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compensator is stabilizing. Because the gradient of the cost function with respect to the 
compensator parameters becomes very large near the stability boundary, the algorithm is 
liwavs forced back towards a stable solution as long as it never steps over the s tab ility 
boundary into an unstable region. 

IV. Two Link robot arm Example 

We used a math model of a two link robot arm (TLA! to demonstrate the capabilities of 
our algorithm. This is the same model used by Yang, and so we were able to verify our 
results by direct comparison. A diagram of the TLA is shown in Figure 4. 

The goal of our design was to control the np position (5) of the arm via a multirate 
compensator. We used the following analog weighting and process noise covariance 
matrices from 1 17]. 


I .21 
0 
0 

- 0 


We assumed perfect measurement and that plant disturbances enter the system 
coincident with the control torques. For the multirate compensators, 9 was sampled and Tj 
was updated every 0.225 seconds: 6 was sampled and T> was updated every 0.028125 
seconds. 

Five different compensators were designed: an analog LQR, a muitirate lead/lead, an 
optimized multirate lead/lead, an optimized multirate general 2 nd order, and an optimized 
single rate general 2 nd order. We used a smooth step input to 6 re f and 0 re j defined as 

follows: 


0 0 0 I 

0 0 0 | 

0 18.5 0 I 

0 0 0 J 


QjJOOl 0 

l 0 0.69444J 


E|ww T J = 


0.69444 0 

0 0.01 


(38) 


Sref(0 — i 




10.01 


t<T c 

t>T c 


®reK0 — 


SrefO) 

Li + L 2 


T c = 0.125 sec 


(39) 


72 



1991 ACC 


Reduced Order Multirate Compensator Synthesis 


11 


and the servo configuration shown in Figure 5 to measure the performance of the different 
compensators [17]. The response of the TLA for the five compensators is shown in 
Figure 6. 

The anaiog LQR compensator used simpie full state feed back. We provided this 
compensator as an example of the response possible using the cost function weighting 
matrices of (38). 

The multirate lead/lead was found using successive loop closures. We designed the 
control loops in the discrete domain so that the eigenvalues of the closed loop system 
matched those we obtained using LQR transformed to discrete time. This compensator 
consists of two simple lead loops: one from 6 to T 2 operating at the fast sampling/update 
rate, and one from 6 to T^ operating at the slow sampiing/update rate. 

The final three compensators were synthesized using our new algorithm and the cost 
.veighung matrices used to design the anaiog LQR compensator. The optimized muidrate 
lead/lead was found by optimizing the pole/zero iocanons and gains of the lead/lead 
compensator found by successive loop closures. 

The optimized muidrate general 2 nd order compensator uses the same sampiing/update 
scheme as the iead/lead compensators but has the compensator structure of (40), where ay, 
bjj, Cy, and dy are the parameters which were optimized. This compensator has the 
maximum number of independent free parameters possible for a second order system. 


0 1 
1 

B J 

1 

b 1 2 1 

c = 

' c n 

Cl2 ' 

D =[ ^ 

di2 

a?2 J 


- t>2l 

1 


. C21 

C22 . 

i d2i 

^22 . 


(40) 


The optimized single rate general 2 nd order compensator is a single rate equivalent of 
the muidrate general 2 nd order compensator. It has the same structure as the muidrate 
general 2 nd order compensator. (40), but uses a single sampling rate which was chosen 
such that the number of computations needed to implement either the muidrate or single rate 
compensators are the same. 

Our results are the same as those obtained using Yang’s algorithm. They demonstrate 
how muidrate compensators can provide better performance than single rate compensators 
by trading lower bandwidth control of the slow modes for higher bandwidth control of the 
fast modes. In this example we were able to reduce the rip response over shoot 40% and 
the peak control torque 25% by using a muidrate controller over a single rate controller. 
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v. Conclusion 

In this paper we have presented an algorithm that can be used to synuiesize reduced 
order muitirate compensators with or without periodically time varying gains. Our 
algorithm is similar to that developed by Yang, except we use an infinite time cost function 
n order to simplify the equations governing the solution. One advantage to this approach 
:s that the expression for the gradient of the cost is relatively simple, making the numerical 
search algorithm more efficient. The disadvantage of this approach is that an initial 
stabilizing compensator is required. However, for most cases a stabilizing compensator 
can be easily obtained using successive ioop closures. 
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APPENDIX 

A. Calculation of the Values of the Lagrange Multiplier 

Given a P k which stabilizes the muitirate system we can calculate the steady 
of Ay. where Ay, is defined by equation (321 rewritten here as (A.l). 

0 = Qi + MP k H + IMP k H) T + (P k H) T Q 2 P k H + Fl k A k+1 F ck - A k 
for k= 1 ,2,...,N with A k - A k+N 


First simplify (A. 11 by defining 


Q 3 =' 


Qi M I 

M T Q 2 j 
I is an identity matrix 


and J k = 1 * 

LP k HJ 


Then (A.l) can be written as 


Ak = JlQsJk + F IkAk+iF c k for k=l,2,...,N with A k = A k+N 

Equation (A.3) represents a periodically time varying Lyapunov equation, 
create an equivalent single rate system by repeated application of (A.3). 

Ak = llkQdJdk + Fj k A k Fdk for k=l,2,...,N with A k = A k+N 


Fdk = F c nc+N-i)Fc(k^-2)Fc(k-*-N-3v • -F C k 


Jdk = [j(k+N-l>F c fk+N-2Fc(k+N-3)- * -F c k | 

J(k+N-2)F c(k-*-N-3)' • -F c k | • • • I Jk] 


Qd = 


Q 3 O .. 

0 03* • • 


0 

0 


0 0 OQ3 


: values 
(A.l) 

(A.2) 

(A.3) 
We can 

(A.4) 

(A.5) 

(A.6) 
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Equation (A.4) is a ume invariant Lyapunov equation which can be solved for A k . 
Once any A k has been found, the propaganon equation (A.3) can be used to find the 
remaining A k . 
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Figure i. A General Multirale Compensator 
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Figure 2. Multirate Sampling Scheme 
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Figure 6a. Tip (5) response to a smooth step input 
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I. INTRODUCTION 

There rue many estaDlished methods for synthesizing muitiratf, compensators (Berg, 
Yang, Mason. Glasson. Amit] but surpnsmgiy few methods for analyzing the robustness 
of these systems. Current robustness analysis methods rely principally on the transfer 
function of the system. A muidrate transfer funcdon. in the traditional sense, does not 
exist, because most muidrate systems are periodically time varying. Without some 
modification, established analysis methods cannot be appiied direcdv to muidrate systems. 

[Thompson! and (Apostoiakis] have both proposed ways to extend existing robustness 
analysis techniques to muidrate systems. Thompson used “Kranc ’ operators to transform a 
special class of muidrate systems, derived from sampled continuous systems, into MIMO 
iingie rate systems. Apostoiakis transformed the general muidrate system into a discrete 
ame single rate MIMO system and then used impulse modulation to produce a continuous 
MIMO system (Boykin's]. In both cases, the inputs and outputs or the new MIMO system 
were comprised of delayed samples of the inputs and outputs of the muidrate system. 
Thompson and Apostoiakis then used muidvariable nyquist criterion and/or unstructured 
singular vaiue analysis to calculate the gain and phase margins for the muidrate system. 

In this paper we will present an alternative approach for extending nyquist criterion and 
singular vaiue analysis to muitirate and periodically time varying systems. Like 
Apostoiakis, we transform the original system into an equivalent time invariant single rate 
system. However, we perform the robustness analysis in the “z” domain. By working in 
the “z” do main we can establish relationships between a muitirate/periodically time varying 
system and its time invariant single rate equivalent. These relationships clarify the 
limitations of nyquist and singular vaiue analysis using single rate equivalent systems. 

The paper is divided into five section. Section I provides some back ground 
information about multirate systems and discusses transfer functions for multirate and 
periodically time varying systems. Section II discusses the application of the nyquist 
stability criterion to these systems; Section III discusses the application of structured and 
unstructured singular values analysis to these systems. Section IV contains a example of 
robustness analysis using structured singular values for a multirate system. Concluding 
re marks follow in Section V. 

I. MULTIRATE AND PERIODICALLY TIME VARYING SYSTEMS 
Before discussing robustness analysis we will first establish the relationship between 
multirate and periodically time varying systems. Then we will define an equivalent single 
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rate system which wtil allow us to combine periodically time varying, muitirate and time 
invariant systems using tradidonai block diagram techniques. 

A general muitirate compensator is shown in Figure 1. Each input (y), output (u), and 
state iT) is sampied/updated at a rate which, in general, represents the desirea bandwidth of 
the input or output with which it is associated, v is the value of v currently available to the 
digital processor from tne zero order hold; while u is the current output from the digital 
processor which is held with a zero order hold to form the output u. A discussion of this 
compensator structure can be found in [Berg 6c Mason, and Yang]. 

Associated with this muitirate compensator is a multirate sampling schedule which 
specifies the samp ling/ update rate for each input, output and state. We define the greatest 
common divisor of ail the sampiing/update periods as the shortest time period (STP) and 
the least common multiple of all the sampling/update periods as the basic time period 
BTP V The integer :V is defined as; 


n = (l) 

STP 

When the sampiing/update rates have ratios which are rational numbers, the 
sampiing/update schedule is periodically time varying and the muldraie compensator can. be 
modeled as a linear periodically time varying system of the form [Mason &c Berg]: 

x(k+l) = A(k)x(k) + B(k)uik) (2) 

yk) = C(k)x(k) + Dk)u(k) (3) 

where Aik) = A(k+N), B(k) = B(k+N), C(k) = C(k+N), and D(k) = D(k+N) 

The sampling period for (2)-(3) is one STP and the period of repeddon is one BTP. 

Any pracdcai muitirate system can be modeled as linear periodically dme varying 
system of the form of (2)-(3). Therefore, we will focus the remainder of the discussion on 
linear periodically time varying systems, of which multirate and single rate are a special 
case. 

Given a periodically rime varying system of the form of (2)-(3), we can create an 
equivalent dme invariant system by repeated application of (2)-(3) over the BTP [Meyer & 
Burrus). The equivalent time invariant system is 

x(N(k+ 1 )) = Ax(Nk) + Bu(Nk) (4) 

y(Nk) = Cx(Nk) + Du(Nk) (5) 

Where 
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A = A(k+N-\)A(k+N-2) 

■ ■ A(k ) 


(6) 

A(N-\)A(N-2)--A(2)B(\)\ A(N-\)A(N-2)---A(2>)B(2) | > 

B(N)] 

(7) 

r a i) 




„ ; C(2)A(1) 

C — l 

! 


(8) 

C(N-\)A(N-2) ■ ■ 

•4(1) J 



r di i) 




C(2)B(1) 




D=i C(3W2)B(1) 




C(N- \)A(N-2)- ■ -A(2)Bn) 

0 

0 

0 


D( 2) 

C( 3)5(2)) 

0 

° i 

(9) 

* 

D(N- 2) 

o ! 


C(N-\)A(N-2)---M2)B(2) 

... C(N-\)B(N-2) 

D(N- 1) j 



where y(Nk) = 

y{Nk) 
y{Nk+ 1 ) 

and u(Nk) = 

u(Nk) 

u(Nk+l) 


\j{Nk+N-\)_ 


_u(Nk+N-\) \ 


Equanons (6) -(10) transforms the linear periodically time varying system. l2)-(3), with 
p inputs, q outputs and a sampling period of one STP to a linear time invariant system, (4)- 
(5), with Np inputs. Nq outputs and a sampling period of one BTP. We will refer to (4)- 
(5) as the equivalent single rate system iESRS) of (2)-(3). It is important to keep in mind 
that the inputs and outputs of an ESRS are comprised of samples of the inputs and outputs 
of a rimff varying system. A consequence of this is the relationship: 

M 2 =imi 2 (11) 

We can calculate the transfer function of (4)-(5) using the following definitions for the 
Z Transform. Let 

oo 

Z(jt(*))=£j:(t)z-* (12) 

i=0 
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Z V UU))=£ .tuN)z*‘ lV 
;=0 

(13) 

then 

y(z iV ,/) = Z N (y{k+l)) 

(14) 

' Zy{y(*)} 

T u v )=1 Zv{y(^+ii} j 
Zy{y(^+N-l)} J 

(15) 

with a similar definitions for mz iV ) and u(z- v ./). 
The transfer funcnon for the £S/?S, (4)-(5), is 


ytz< v ) = Gv(z ,v '):t(z' v ) 

where Gy(z ,v ) = Cl/z^ - A)' l B + D 

(16) 

(17) 

The transfer funcdon is written as G/v(z lV ) to emphasize that the sampling period of the 
ESRS is one BTP . or (V times the sampling period of the time varying system t2)-(3). 

So far. the ESRS has only been applied to periodically time varying systems. We 
could, however, calculate the ESRS of a time invariant system - in this case N can be any 


integer. 

The transfer function for the ESRS of a time invariant system can be calculate using 
(17). Alternatively, Gy(z ,v ) can be calculated directly in terms of the transfer function of 

the time invariant system, G(z). Given 

y(z) = G(z)u(z) (18) 

and following [Meyers & Burras) we can write 


y(z) = y(z N ,l) 
1=0 ^ 

™ i=0 

irv 

where <p =e * 

Combine ( 18)-(20) to obtain 


(19) 

( 20 ) 
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V-l AM 

v(z' v ,/) = z* z'" 1 # " m G(z0‘)u(z ,s, .m) 

y 1=0 m=u 


( 21 ) 


From 1 15) ana (21) the F 71 row and m ;,<l column or Gk/(z'') is given by 


v-i 


Ul-m) , 


= "CMl 

iV :=0 


( 22 ) 


For a ume invanant system. G’v(z ,v ) is made up of time and frequency shifted versions 
of G(z). A special case of (22) occurs when G(z0 l ) = G(z) for i - 0,1,...„V-1 . 


If G(z0) = G(z) for t=0,l,...„V-l then G N (z^) itn = 0 ll cr ^ c (23) 

The simDiest G(z) satisfying (23) is G(z) = constant. Equanon (23) is an important 
reiauonsniD wnicn wiii be used in Secnons II and HI. 

Equauons (6)-(10), (16)-( 17) or (21) can be used to compute state space and transfer 
function aescnptions for the ESRS of a periodically time varying or time invanant system. 
The aavantage of the ESRS is that is it allows us to manipulate dme invariant and 
periodically time varying systems (e.g. muitiratc) as if they were both time invariant. The 
state space or transfer functions descriptions can be used to calculate input-output relations 
for systems in series or in a feedback loop just as in classical control [Khargonekar]. In 
addition. [Konoj has shown that if the ESRS is stable then the time varying system from 
which it was derived will be stable. So, we need only worry about the stability of the 
ESRS. 


II. Nyquist Stability Criterion 

We can determine the stability of the periodically time varying system in Figure 2 by 
applying standard multiloop nyquist [McFarlane ...] criterion to the ESRS, since the 
periodically time varying system will be stable if its ESRS is stable. The ESRS return 
difference is given by 


/ - G N (z")A N {z") ( 24) 

and the nyquist contour is 

z ,v - gjot 0£ w £ 2ic 

When the periodically time varying system is SISO, we can determine traditional gain 
and phase margins from the nyquist plot. Recall that when A(z0‘) = A(z), A^(z^) is a 
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diagonal transfer function matrix wuh A(z) on the diagonal. Thus, tne ESRS for Figure 2 
with gain and phase uncertainty can be written as 

Gy(z lV 

^actual — Gy(Z‘’ ) nominal (25) 

Phase ana gam margins from me nyquist plot can be interpreted in the traditional sense 
even tnougn the ESRS is MIMO because the inputs and outputs are correlated in time and a 
constant earn applies equally over ail time. [Thompson] amveti at mis same results using 
Franc operators. 

When me ume varying system is MIMO. the standard MIMO nyquist restrictions apply. 
For MIMO time varying systems, it is best to use a norm based approach sucn as singular 
value anaivsis. 


in. Structured and Unstructured Singular Values analysis 
In the previous section we saw that the muitiloop nyquist stability criterion can be 
appiied to an ESRS to determine the stability of a periodically time varying system. In this 
section we will see that, with some limitations, both structured and unstructured singular 
value analysis can be appiied to the ESRS to determine the robustness properties of a 
periodically time varying system. 

Given stable transfer functions G(z) and A(z) it has been shown that the closed loop 
system will remain stable through out continuous changes in A(z) if 

del (/ - G(z)A(z)) *0 (26) 

or g (/ - G(z)A(z)) > 0 (27) 

is satisfied around the nyquist contour, subject to certain restriction on Gz) and A(z) 
[Maciejousky,...]. By direct application of (26)-(27), a periodically time varying system 
will be stable if 


det U - Gv(z‘ v )A v(z' v )) * 0 (28) 

or g (/ - Gy(z ,v )Ay(z' v )) > 0 (29) 

is satisfied around the nyquist contour, because a periodically time varying system will be 
stable if its ESRS is stable. From (28)-(29) it follows that most singular value robustness 
tests can be appiied directly to a ESRS to determine the robustness properties of a 
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periodically time varying system. The resuits, though, must be interpreted in light of the 
fact mat some or the inputs and outputs of the ESRS are time correlated. In the following 
r aragrapns we wiil discuss the imoortant differences between a single rate system and an 
ESRS ana how tnese affect singular value analysis. 

1) There are qN, not q, singular vaiues associated with each point on the nvquist 
sontour tor tne ESRS of a time varying system with oniy q inputs and ourputs. The 
jadirionai singular vaiues come from the time correlated inputs and outputs of the ESRS. 
Rememoer that the sampling penod of the ESRS is one BTP. N dmes slower than the rime 
varying system from which it was derived: but the ESRS has N rimes as many inputs and 
outputs as the original rime varying system. The key point is that ail of these singular 
values are important in determining the roousmess of a periodically time varying system. 

if an ESRS is generatea from a time invariant system, the singular values of the ESRS 
and the singular values of the onginai single rate system are related by the following 
expression. 


aGrteJN*) - 1 oG(*°e>"), oG(tf> ■ ■ • <sG(<t> NA e^)\ (30) 

In (30), singular vaiues associated with frequencies above 1 /BTP in G(z) are reflected 
back to lower frequencies in Gy(z iV ). It follows directly from ( 30) that 

l|C K ( 2 'V)|L-!|G(z)|_ ,31 > 


where j|G(z)L=sup a[G(£/ w )j 

2) The ESRS imposes a structure on any uncertainly Ay. Any A v in (28)-(29) must 
obey (17) or (22); this automatically imposes a structure on A v .. The problem is, we are 
often interested in a time invariant plant uncertainty. A, that destabilizes the system and not 
in A v . a (Ay) found using unstructured singular value analysis is often overly conservative 

because it accounts for not oniy the fictitious perturbations normally associated 
unstructured singular values but also for rime varying and non-causai perturbations. It is 
important to remember that o( Ay) found using unstructured singular values can be 
extremely conservative and may not reflex 5( A). 
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Structurea singular values provides a mechanism tor rinding A. For an ESRS with a 
time invariant uncertainty. A(z), the definition of the structured singular vaiue. y., can be 

written as 


i 0 if detii - Gv(z ,v )Av(z' v )) * 0 for any A e A B d 
j. i q s _ i-i (32) 

v ' min tC(A(z)): det (7 - Gv(z)Av(z)) =0 omerwise 

\ Ae abo » 

Ann is the form of the permissible block diagonal perturbations A: and the structure of Ay 

lj 

must satisfy equation (22). 

Unfortunately tor a general A(z), the structure of A^(z) is otten very complex and 
finding a good estimate of size or A is difficult. However, wnen A is a constant, as is the 
case for many proolems. 

A v = diag( A, A. ... A) with N blocks (33) 

and structured singular vaiue analysis can be used to determine A. 

When A is a nme varying, but not a function of “z”, A w becomes 

A^ = diag( A(l), A(2), ... A(N)) (34) 

and has no repeated blocks. Equation (34) must be interpreted with care - (34) implies that 
the vaiue of Aik) is constant over the sampling interval, STP, and changes instantaneously 
to Aik+l) at the next sampling instant. This may not be a good model of time varying 

uncertainty. 

3) When A is a constant then each A block of A N can be scaled independently. Using 
the block diagonal scaling property of |i [Maciejowski], and (33) or (34) it is 
straightforward to see that 

[f A is a p by p matrix then \i{DGy(z‘V)D‘ l ) = | 1 (Gn(z' v )) (35) 

where D = id x l p> dJ p , ... d^ p ) and / is a p by p identity matrix 

In addition if A is block diagonal then each of sub-block of A can be scaled in a similar 
manner. 

An interesting result of (35) is that the upper bound for y.(G/v(z‘ v )) given by 
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i uG/v(z- v )) < i^DGNiz^D' 1 ) (36) 

:s the same whether A is periodically time varying or time invariant subject to the 
interpretation of a time varying A mentioned in item 2. The upper bound for the time 
invariant case found using (36) is of course more conservative. 

4} Singular vaiue piots of ESRS transfer function matrices should not be interpreted in 
:he frequency domain. The ESRS has time correlated inputs and outputs. The response 
from one input to one output represents only pan of the total signal between the input and 
output of the penodicallv time varying system. 

A meaningful auanuty for an ESRS is its infinity norm. From (11) and (Francis) we 

jon write that 


sup = sup = j|Gy(z /V )|| „ for i|u|| < oo (37) 

ml 

Thus, ilG*(z*)L can be interpreted as the maximum gain of the system for all u with a 

bounded two norm, just as in the single rate case. For the single rate case the maximum 
gain occurs when u is sinusoidal - this is not necessarily true for the penodically time 
varying system. 

Singular vaiue analysis of an ESRS. both structured and unstructured, can be used to 
determine the robustness of periodically time varying system. As we have discussed, there 
are limitations to this analysis because the inputs and outputs of an ESRS are time 
correla t ed. 


IV. Two link arm Example 

The results of the previous sections will be illustrated by calculating the gain margins 
for a planar two link robot arm (TLA) using structured singular values. Two different cases 
are considered: 1) the TLA with a 2 nd order muitirate compensator and 2) the TLA with a 
2 nd order single rate compensator. 

The TLA is shown in Figure 3 and is described further in (Berg and Yang). The two 
compensators were designed to minimize a cost function quadratic in the states and controls 
us ing the optimization method described in [Mason & Berg]. 
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The 2 na order muitirate compensator uses the compensator structure ot (38), where fl.y, 
:,.. f c,„ and d ;; are the parameters which were optimized. 

I ■ Hi 0 ~ ! b=i 1 bl2 'l C = t C 1 1 C!2 1 D=i dn dl - 2 'i (38) 

' _ 0 U 2 2 J . b 2 l 1 ] - C21 022 J - d 21 d 22 1 

The sampiing/update rates tor the compensator are listed in Table 1. In addition the 
compensator state associated with 8 and Tj is updated at the slow rate while the state 
associated with 6 and To is updated at the fast rate. For the muitirate compensator STP — 

.028125. BTP = 225 andiV = 8 


i 

i 

Multirate Comoensator i 

Single Rate Compensator 

9 

0.225 s 

0.05 s 

5 ! 

0.028125 s 1 

0.05 s 

i 

T ! 

M ! 

0.225 s 

0.05 s 

1 i 

| 0.028125 s 

0.05 s 


Table 1. Sampiing/Update Periods for the Compensators 


The single rate 2 nd order compensator is the single rate equivalent of the muitirate 
compensator. It has the same structure as the 2 nd order multirate compensator and 
minimizes the same cost function, but uses a single sampling/update rate. This 
samp iin g/update rate was chosen such that the number of computations required to 
implement either the muitirate or single rate compensators during real-time operation is the 
same. The sampiing/update rate for the compensator is shown in Table 1 . 

A block diagram of the TLA. compensator, and output gain uncertainty is shown in 
Figure 4. The bock diagram in Figure 4 can be cast into the standard structured uncertainty 
model shown in Figure 5 where the gains and hy are allowed to vary independently. An 

upper bound on the structured singular values for the multirate and the single rate cases was 
calculated using the following [Safonov). 

q(Q) < irtfc(DQP A ) < k P (Q) (39) 

where k p (Q) is the Perron-Frobenius eigenvalue of Q 
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For the muidrate case, the ESRS for the piant. compensator and uncertainty was 
calculated for N = 8 . They were combined as shown in Figure 5 and an upper bound for 
a, as z traversed the nvquist contour, was calculated using (39). A lower bound on the 

maximum singular value of the gain matrix is given in Table 2. 

For the sineie rate case an upper bound on p. was calculated using two different 
methods. First, \i was calculated direcdv for the single rate system. N = l . An exact 
va i ue for u can oe calculated because there are only two blocks in the uncertainty matrix 
Doyiej. Next, an ESRS was constructed for the single rate case using N = 8. An upper 
bound for u was calculated using (39). For the ESRS system, the uncertainty matrix has 8. 

2 by 2 blocks. These resuits are summarized in Table 2. 


Design 

, . ffc, 0 1 11 

Gam Margin 01 q ^ ( = — 

Muidrate 2 nd Order 

0.535 j 

Single Rate 2 nd Order 

! 0.513 

i 

Sineie Rate 2 nd Order usine ESRS 

0.389 ' 


Table 2. Gain Margins for TLA with Multirate and Single Rate Compensator 


The two p. es timat es for the single rate case illustrate the disadvantage of using (39) to 
calculate the upper bound of p. for an ESRS systems. As in the muidrate case, the ESRS 
single rate case accounts for periodically time varying uncertaindes. resulting in a 
conservative estimate of u. See item 3, Section IIL 

For the ass umed uncertainty model the muidrate compensator was siighdy more robust, 
even given the conservativeness of the estimate for (L The muidrate compensator is able to 
compensate for larger gain uncertainty because it has higher bandwidth control of the 
second link than does the single rate compensator. 


v. Summary and Conclusions 

In this paper we have shown how nyquist criterion and singular value analysis can be 
applied to muidrate and periodically time varying systems using their ESRS. For SISO 
systems, traditional gain and phase margins can be found by direct application of the 
nyquist criterion to the ESRS. For MIMO systems, structured singular values can be used 
to determine the maximum size of an uncertainty. The results of singular value analysis, 
though, must be interpreted in light of the fact that some of the inputs and outputs of the 
ESRS are time correlated. We pointed out several important resulting l imit a t i o ns of 
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mauiar vaiue analysis using an ESRS. Finally we demonstrated robustness analysis for a 
wo iink arm wun a muiurate compensator usmg structured singular value. 
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APPENDIX 

Lemma: oGy(eJ Nw ) = i l oG(0°c>), oG(0 l eJ a) ), ■ ■ • aG(0 NA ei Qi )\ 

Proof: From the definition of a transfer function we can write 

y(z) = G(z)u(z) 


(A.l) 


where v(z) = 


" ,,0.1 
y{0 z) 


y(0 l z) 

! 

, w(z) = 

, ,'V-l % 

iy(0 -). 



U(0°Z) j 
U(0 1 Z) | 


u(0 N ' l :)j 


G(z) =1 


G(0 z) 0 
0 G(0 1 z) 


0 


V-l 

0 0 0 G(0 z) J 


(A.2) 


N - 1 


From 1 12), y(z) = ^ z ml y(z N ,l) so that we can write y(z) = Ty(z N ) 

1=0 


where T = 


/ z-l/ ... z*<A/-l>/ 

/ (<f> l z)- 1 / 


i 


(A.3) 


|_ / (^z)' 1 / ... (0 iWA z)W)[ J 

T has the property that IT* =NI if z is evaluated on the unit circle and / is an identity 
matrix of appropriate dimensions. 

Then 

y(z N ) = r 1 G(z)Tu(z") so that G w (z*) = r l G(z)T (A.4) 


Now using the fact that a 2 (A) equals the eigenvalues of A* A, and that the eigenvalues 
of a block diagonal matrix are the eigenvalues of the individual block it follows that 


= [<jG(0 O eJ a \ aG(0 ■ . . aG(0 S ' 1 eJ°>)} (A.5) 
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:(&+l) = Az(k) + By(k) 
u(k) = Cz(k) + Uy(k) 


ZOH I- 


ZOH t- 


ZOH H 


Digital Processor 


Figure 1 . A General Multirate Compensator 


A(r) 


- G(t) 


Figure 2. Periodically Time Varying System 
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Parameters: Mass Length 

L, 1.235 kg 0.965 m 

L 2 0.163 kg 0.167 m 

Inputs: Torque Tj and T 2 
Outputs: 6 and 6 


Figure 3: Diagram oi the Two Link Arm 


TLA Plant 


k,j 

Gain Uncertainty 




Figure 5. Structured Uncertainty Model for TLA 
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